% number of elements in table
t_numel = num_unit*num_drug*2;
% compute attention averages
t_rate(:,:,1) = mean(rate_summary.rate(:,:,idx_attention_cond==1),3);
t_rate(:,:,2) = mean(rate_summary.rate(:,:,idx_attention_cond==0),3);
t_FF(:,:,1) = mean(rate_summary.FF(:,:,idx_attention_cond==1),3);
t_FF(:,:,2) = mean(rate_summary.FF(:,:,idx_attention_cond==0),3);
t_PSTH(:,:,1,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==1,:),3);
t_PSTH(:,:,2,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==0,:),3);
t_gain = arrayfun(@(x)(x.rhat),rate_ROC.gv,'UniformOutput',true);
t_gain_log = log(t_gain);
% use idx_attention_roc for ROC and MI, not that this is different from the
t_roc_attend(:,:,1) = mean(rate_ROC.roc_attend(:,:,idx_attention_roc==1),3);
t_roc_attend(:,:,2) = mean(rate_ROC.roc_attend(:,:,idx_attention_roc==0),3);
t_mi_attend(:,:,1) = mean(rate_ROC.mi_attend(:,:,idx_attention_roc==1),3);
t_mi_attend(:,:,2) = mean(rate_ROC.mi_attend(:,:,idx_attention_roc==0),3);
% create metadata columns
[t_att] = deal(zeros(size(t_rate)));
t_unit = categorical((1:num_unit)', 'Ordinal', false);
t_unit = repmat(t_unit, [1, num_drug, 2]);
t_drug_on = repmat(t_drug_on, [num_unit, 1, 2]);
t_unit = reshape(t_unit, [t_numel, 1]);
t_drug_on = reshape(t_drug_on, [t_numel, 1]);
t_att = reshape(t_att, [t_numel, 1]);
[t_att_label,t_drug_on_label] = deal(cell(length(t_att),1));
t_att_label(t_att==1) = {'Attend RF'};
t_att_label(t_att==0) = {'Attend away'};
t_drug_on_label(t_drug_on==0) = {'Drug off'};
t_drug_on_label(t_drug_on==1) = {'Drug on'};
% reshape matrices to single columns
t_rate = reshape(t_rate, [t_numel, 1]);
t_FF = reshape(t_FF, [t_numel, 1]);
t_roc_attend = reshape(t_roc_attend, [t_numel, 1]);
t_mi_attend = reshape(t_mi_attend, [t_numel, 1]);
t_gain = reshape(t_gain, [t_numel, 1]);
t_gain_log = reshape(t_gain_log, [t_numel, 1]);
t_gv = reshape(t_gv, [t_numel, 1]);
t_PSTH = reshape(t_PSTH, [t_numel, length(rate_PSTH.time)]);
t_drug_on_label, t_att_label, ...
t_rate, t_FF, ... % rate, FF
t_roc_attend, t_mi_attend, t_gain, t_gain_log, t_gv, ... % roc, MI, gain, gain_log, GainVariability
'VariableNames', {'unit', 'drug_on', 'attention', 'rate', 'FF', 'roc', 'MI', 'gain', 'gain_log', 'gv', 'PSTH'});
att_table = join(att_table, unitlist, 'Keys', 'unit');
% define categorical variables
att_table.drug_on = categorical(att_table.drug_on, 'Ordinal', false);
att_table.attention = categorical(att_table.attention, 'Ordinal', false);
% reorder categories to label_* order
att_table.drug_on = reordercats(att_table.drug_on, label_drug_onoff);
att_table.attention = reordercats(att_table.attention, label_attention);
% define combined classes
att_table.attention_drug_on = att_table.attention .* att_table.drug_on;
att_table.attention_drug_on = categorical(regexprep(cellstr(att_table.attention_drug_on), ' Drug', ', Drug')); % add a comma in between classes
% clear temporary variables
disp(head(att_table))
unit drug_on attention rate FF roc MI gain gain_log gv PSTH Subject Task Date Drug EjectCurrent Weight peak_to_trough_time unit_class waveform EjectCurrent_centered s_stim s_att s_dru s_dir s_att*dru s_att*dir s_dru*dir s_att*dru*dir attention_drug_on
____ ________ _________ ______ ______ _______ ________ ________ ________ _____________________ _______________ _______ _____ __________ ________ ____________ ______ ___________________ __________ ______________ _____________________ __________ __________ __________ __________ _________ _________ _________ _____________ ___________________
1 Drug off Attend RF 8.2778 2.1547 0.6376 0.12527 0.14343 -1.9419 [1×1 GainVariability] [1×1251 double] W gratc 2015-05-21 Dopamine 25 9.29 345.6 Broad [1×187 double] -34.25 0.6978 0.0023999 0.57272 0.83275 0.55229 0.52515 0.30834 0.64143 Attend RF, Drug off
2 Drug off Attend RF 8.8506 1.4285 0.70263 0.17977 NaN NaN [1×1 GainVariability] [1×1251 double] W gratc 2015-05-21 Dopamine 25 9.29 205.2 Narrow [1×187 double] -34.25 0.79882 5.7922e-14 0.16206 0.25224 0.027989 0.54777 0.21153 0.20304 Attend RF, Drug off
3 Drug off Attend RF 23.394 3.7531 0.66866 0.12168 0.11193 -2.1899 [1×1 GainVariability] [1×1251 double] W gratc 2015-05-21 Dopamine 25 9.29 388.8 Broad [1×187 double] -34.25 0.05165 4.2769e-08 0.45615 0.31594 0.054758 0.42777 0.44512 0.26137 Attend RF, Drug off
4 Drug off Attend RF 11.453 1.8759 0.64466 0.11753 0.065949 -2.7189 [1×1 GainVariability] [1×1251 double] W gratc 2015-05-28 Dopamine 60 NaN 394.2 Broad [1×187 double] 0.75 0.29917 1.5365e-08 0.98854 0.90503 0.1141 0.97673 0.61297 0.38682 Attend RF, Drug off
5 Drug off Attend RF 28.963 5.2381 0.61176 0.088129 0.18975 -1.662 [1×1 GainVariability] [1×1251 double] W gratc 2015-05-29 Dopamine 80 NaN 183.6 Narrow [1×187 double] 20.75 2.0498e-22 6.6763e-07 0.0087419 0.88981 0.43338 0.71566 0.67429 0.44351 Attend RF, Drug off
6 Drug off Attend RF 44.259 7.4757 0.73118 0.19062 0.13876 -1.975 [1×1 GainVariability] [1×1251 double] W gratc 2015-06-03 Dopamine 75 9 302.4 Broad [1×187 double] 15.75 5.3262e-45 2.8411e-14 0.00011895 0.74406 0.24956 0.68049 0.2666 0.1016 Attend RF, Drug off
7 Drug off Attend RF 24.173 2.6038 0.73431 0.1636 0.065311 -2.7286 [1×1 GainVariability] [1×1251 double] W gratc 2015-06-03 Dopamine 75 9 329.4 Broad [1×187 double] 15.75 1.8937e-29 2.2529e-13 0.011994 0.00099457 0.80152 0.022988 0.93706 0.83731 Attend RF, Drug off
8 Drug off Attend RF 12.152 5.0028 0.57967 0.14417 0.33173 -1.1034 [1×1 GainVariability] [1×1251 double] W gratc 2015-06-04 Dopamine 90 8 291.6 Broad [1×187 double] 30.75 7.4239e-14 0.00078848 0.075078 0.043667 0.91433 0.34239 0.58849 0.20115 Attend RF, Drug off
% save to csv file for R analysis
table_selectivity = get_unit_selectivity(att_table, selectivity_criterium);
R_table = att_table(table_selectivity,:);
R_table = removevars(R_table,{'PSTH','waveform'});
writetable(R_table, fullfile(path_population,sprintf('population_attention_%s.csv',selectivity_criterium)))
[colors, labels] = get_colors('att_drug');
colors = reshape(colors, [4,3]);
labels = reshape(labels, [4,1]);
linestyle = {'-','-','-','-'};
events = {'stim', 'cue', 'dim'};
label_event = {'Stimulus onset', 'Cue onset', 'First-dimming'};
xlabel_event = {{'Time from','stimulus onset (ms)'}, {'Time from','cue onset (ms)'}, {'Time from','first-dimming (ms)'}};
xlim2use = {[-200 750], [-100 750], [-1000 200]};
ncol = length(label_event) + 2;
nrow = length(label_drug);
idx_subplot = {1,2,3,6,7,8,[4 5 9 10]};
[fH, fSet] = plotj_initFig('width', 26, 'height', 11, 'Journal',plot_conventions);
fSet.subplotGap2 = fSet.subplotGap .* [1.3 0.30];
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(irow)});
for ievent = 1:length(label_event)
rate_PSTH = get_population_data(recordinglist, 'spike_rate_PSTH', path_data, events{ievent});
% convert samples to att_table format
t_numel = num_unit*num_drug*2;
t_PSTH(:,:,1,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==1,:),3);
t_PSTH(:,:,2,:) = mean(rate_PSTH.samples(:,:,idx_attention_cond==0,:),3);
t_PSTH = reshape(t_PSTH, [t_numel, length(times)]);
subtightplot(nrow, ncol, idx_subplot{iplot}, fSet.subplotGap2, fSet.subplotMargin, fSet.subplotMargin)
plotj_initAx(fSet, 'axlabel', irow, 'axlabelDisplacement', [0.07, 0.00]);
if ievent==(ceil(length(label_event)/2))
ht = title([label_drug{irow} label_drug_ext{irow}], 'FontSize', fSet.Fontsize_title);
% ht.Position = ht.Position + [0.05 0 0];
axh.YAxis.Visible = 'off'; % remove y-axis
% axes('Color','none','YColor','none')
% plot each attention and drug condition
conditions = categories(att_table.attention_drug_on);
for iatt_drug = 1:length(conditions)
idx_att_drug = att_table.attention_drug_on==labels{iatt_drug};
idx = idx_unit & idx_att_drug;
h(icol) = boundedline(times, ...
squeeze(mean(t_PSTH(idx,:),1)), ...
squeeze(std(t_PSTH(idx,:),[],1))/sqrt(length(find(idx))), ...
'cmap', colors(iatt_drug,:), 'alpha');
set(h(icol), 'LineWidth', fSet.LineWidth)
set(h(icol), 'LineStyle', linestyle{iatt_drug})
ylabel('Normalized firing rate (%)', 'FontSize', fSet.Fontsize_text)
text(200, 30, sprintf('n=%d', length(find(idx))))
hleg = legend(h, labels, 'FontSize', fSet.Fontsize_text, 'Box', 'Off');
hleg.Position = hleg.Position + [0.05 0.15 0 0];
xlabel(xlabel_event{ievent}, 'FontSize', fSet.Fontsize_text)
hleg = legend(h, labels, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'Location', 'NorthWest');
hleg.Position = hleg.Position + [-0.04 0.53 0 0];
% plot example cell, plot raster for example cell, only drug/attention
axH = subtightplot(nrow, ncol, idx_subplot{iplot}, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', 2, 'axlabelDisplacement', [0.07, 0.00]);
[colors] = get_colors('att_drug');
% load data for single unit
rec2plot = {'W','2015-10-01'};
idx_unit = unitlist.Subject==rec2plot{1} ...
& unitlist.Date==rec2plot{2} ...
& unitlist.unit==categorical(unit2plot(1));
loadfilename = fullfile(path_data, string(unitlist.Subject(idx_unit)), string(unitlist.Date(idx_unit)), 'unit.mat');
unit = load(loadfilename);
loadfilename = fullfile(path_data, string(unitlist.Subject(idx_unit)), string(unitlist.Date(idx_unit)), 'trialdata.mat');
% [trialdata, unit] = remove_excluded_trials(trialdata, unit, unit2plot);
% select only trial window for which this unit has spikes
% raster for each condition
idx_attend = idx_attention_cond([trialdata.cond_num])';
idx_drug = [trialdata.drug]' + 1;
for iatt = 1:length(unique(idx_attention_cond))
& abs(idx_attend-2)==iatt;
[xspikes, yspikes] = spike_raster( sps(:,trial_index), xlim2use{event2plot}, unit2plot(2));
y = squeeze(mean(rate_PSTH.samples(idx_unit,idrug,abs(idx_attention_cond-2)==iatt,:),3));
h(icol) = plot(x, y, 'Color', colors(iatt,idrug,:), 'linew', fSet.LineWidth, 'LineStyle', linestyle{icond});
plot(xspikes, yspikes+yoffset+ymax, 'Color', colors(iatt,idrug,:), 'linew', 1.5)
% set y-offset for next iteration
ymax = ymax + max(yspikes(:));
plot([-500 0], [10 10], 'Color', [0.5 0.5 0.5], 'linew', 5)
set(gca,'YTick',[0 yoffset],'YTicklabel',[0 round(rate_PSTH.maxhist(idx_unit)*1000)]);
xlabel('Time from first-dimming (ms)', 'FontSize', fSet.Fontsize_text)
ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)
ylim2use = yoffset+ymax+80;
p = [unitlist.s_att(idx_unit) unitlist.s_dru(idx_unit) unitlist.("s_att*dru")(idx_unit)];
[pstring,starstring] = get_significance_strings(p, 'rounded', 0, 'factorstring', {'attention', 'drug', 'interaction'});
text(axH, -750, ylim2use*0.95, pstring, 'FontSize', fSet.Fontsize_text)
data2plot = {'rate','FF','gain_log'};
colors = get_colors('spikewidth');
[colors_violin, labels_violin] = get_colors('att_drug');
colors_violin = reshape(colors_violin, [4,3]);
labels_violin = reshape(labels_violin, [4,1]);
ncol = length(data2plot) + 1 + 1; %
[fH, fSet] = plotj_initFig('figNum', 2, 'width', 26, 'height', 13, 'Journal',plot_conventions);
fSet.subplotGap = [0.1, 0.06];
fSet.subplotMargin = fSet.subplotMargin / 2;
fprintf('Testing model: %s \n', label_drug{idrug})
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
num_unit_plot = length(unique(att_table.unit(idx_unit)));
fprintf('\t Found %d units \n', num_unit_plot)
for idp = 1:length(data2plot)
subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);
axH = plotj_initAx(fSet);
% v_att_table = unstack(att_table(idx_unit,:), data2plot{idp}, {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {data2plot{idp},'drug_on','attention','unit','attention_drug_on'}));
case {'gain_log', 'gain'}
gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit'}, 'VariableNamingRule', 'preserve');
variable_names = gv.Properties.VariableNames;
for iatt_drug = 1:length(labels_violin)
mus = [gv.(labels_violin{iatt_drug}).mus];
s2s = [gv.(labels_violin{iatt_drug}).s2s];
rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);
gv_pop = GainVariability;
m2v = gv_pop.pred_variance(xval_gain);
scatter(mus, s2s, fSet.MarkerSize, colors_violin(iatt_drug,:), 'filled', 'MarkerFaceAlpha', 0.5);
plot(m2v.mean, m2v.variance,'Color', colors_violin(iatt_drug,:), 'linew', fSet.LineWidth);
xlim([0.8 max(xval_gain)])
ylim([0.8 max(xval_gain)])
plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)
set(gca,'xscale','log','yscale','log')
% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)
if idrug==length(label_drug)
xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)
ylabel('Variance (spikes^{2})','FontSize',fSet.Fontsize_text)
subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);
axH = plotj_initAx(fSet);
v = violinplot(v_att_table.(data2plot{idp}), v_att_table.attention_drug_on, 'ShowMean', true);
XTickLabel = get(gca,'XTickLabel');
v(iv).ScatterPlot.SizeData = 10;
idx = strcmpi(XTickLabel(iv), labels_violin);
v(iv).ViolinColor = colors_violin(idx,:);
% plot lines for each unit
[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));
% data_x(:,iv) = v(iv).ScatterPlot.XData;
% data_y(:,iv) = v(iv).ScatterPlot.YData;
if ~any(isnan(v_att_table.(data2plot{idp})))
% works if data is complete
data_x(:,iv) = v(iv).ScatterPlot.XData;
data_y(:,iv) = v(iv).ScatterPlot.YData;
% catch if nan values are present
idx = strcmpi(XTickLabel(iv), labels_violin);
idx_trial = v_att_table.attention_drug_on==labels_violin{idx};
data_y(:,iv) = v_att_table.(data2plot{idp})(idx_trial);
plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);
lme_population = fitlme(v_att_table,[data2plot{idp} ' ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
coefficients = lme_population.Coefficients;
betas = coefficients.Estimate(2:end);
betas_se = coefficients.SE(2:end);
p_val = coefficients.pValue(2:end);
% p_string = get_significance_strings(p_val, 'rounded', 0);
p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});
p_string{end+1} = sprintf('n=%d',num_unit_plot);
% % full betastring: Estimate, SE, p-value
% sprintf('drug: %s = %1.2f%s%1.2f, %s', '\beta', betas(1), '\pm', betas_se(1), p_string{1}), ...
% sprintf('att: %s = %1.2f%s%1.2f, %s', '\beta', betas(2), '\pm', betas_se(2), p_string{2}), ...
% sprintf('att %s drug: %s = %1.2f%s%1.2f, %s', '\times', '\beta', betas(3), '\pm', betas_se(3), p_string{3}), ...
% sprintf('n = %d',length(unique(v_att_table.unit)))};
% % part betastring: Estimate, SE, p-value
% sprintf('drug: %s = %1.2f, %s', '\beta', betas(1), p_string{1}), ...
% sprintf('att: %s = %1.2f, %s', '\beta', betas(2), p_string{2}), ...
% sprintf('att %s drug: %s = %1.2f, %s', '\times', '\beta', betas(3), p_string{3}), ...
% sprintf('n = %d',length(unique(v_att_table.unit)))};
case {'gain_log', 'gain'}
set(gca, 'ylim', YLIM .* [1 1.5])
set(gca, 'ylim', YLIM .* [1 1.15])
x_pos = get_value_range(tmp_x, 0.95);
y_pos = get_value_range(tmp_y, 0.95);
h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right', 'FontSize', fSet.Fontsize_text);
p_masked = [p_masked ; false];
plotj_text_emphasise(h_text, 1, 'italic', p_masked);
plotj_text_emphasise(h_text, 1, 'bold', p_masked);
ylabel({label_drug{idrug}, 'Firing rate (Hz)'}, 'FontSize', fSet.Fontsize_text)
ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)
ylabel('Fano Factor', 'FontSize', fSet.Fontsize_text)
ylabel('\sigma_{G}^{2}', 'FontSize', fSet.Fontsize_text)
ylabel('log(\sigma_{G}^{2})', 'FontSize', fSet.Fontsize_text)
set(gca,'XTickLabelRotation',25)
set(gca, 'XTickLabel', '')
subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabel', iplot + plotlabel_offset, 'axlabelDisplacement', [0.05, 0.02]);
axH = plotj_initAx(fSet);
idx_cond = att_table.attention=='Attend RF';
idx = idx_unit & idx_cond;
tmp_data = zeros(length(find(idx))/2,1);
tmp_data(:,1) = att_table.roc(idx & att_table.drug_on=='Drug off');
tmp_data(:,2) = att_table.roc(idx & att_table.drug_on=='Drug on');
plotj_scatter(tmp_data, ...
'markerStyle', markerstyle, 'MarkerSize', markersize, ...
'markerFaceColor', [0.5 0.5 0.5], 'markerFaceAlpha', 0.5, ...
'markerEdgeColor', [0.5 0.5 0.5], 'markerEdgeAlpha', 0.5, ...
P(idrug) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');
p_string = get_significance_strings(P(idrug), 'rounded', 0);
delta_data = tmp_data(:,2)-tmp_data(:,1);
d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');
fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
std(delta_data)/sqrt(size(tmp_data,1)), ...
x_pos = get_value_range(tmp_x, 0.05);
y_pos = get_value_range(tmp_y, 0.95);
h_text(idrug) = text(x_pos, y_pos, ...
sprintf('%s (n=%d)', p_string, size(tmp_data,1)), ...
'FontSize', fSet.Fontsize_text);
xlabel('Attention AUROC no drug','FontSize', fSet.Fontsize_text)
ylabel('Attention AUROC drug','FontSize', fSet.Fontsize_text)
end
Linear mixed-effects model fit by ML
Model information:
Number of observations 124
Fixed effects coefficients 4
Random effects coefficients 31
Covariance parameters 2
Formula:
rate ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
848.12 865.04 -418.06 836.12
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 19.132 2.7502 6.9564 120 1.9744e-10 13.687 24.577
{'drug_on_Drug on' } -2.3639 0.38775 -6.0964 120 1.3585e-08 -3.1316 -1.5962
{'attention_Attend RF' } 2.6939 0.38775 6.9476 120 2.0643e-10 1.9262 3.4616
{'drug_on_Drug on:attention_Attend RF'} -0.2884 0.38775 -0.74379 120 0.45846 -1.0561 0.47931
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 15.16 11.76 19.543
Group: Error
Name Estimate Lower Upper
{'Res Std'} 4.3178 3.7398 4.9851
Linear mixed-effects model fit by ML
Model information:
Number of observations 124
Fixed effects coefficients 4
Random effects coefficients 31
Covariance parameters 2
Formula:
FF ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
447.27 464.19 -217.63 435.27
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 3.4901 0.21273 16.406 120 1.9163e-32 3.0689 3.9112
{'drug_on_Drug on' } -0.02305 0.10547 -0.21855 120 0.82737 -0.23186 0.18576
{'attention_Attend RF' } 0.10797 0.10547 1.0238 120 0.30801 -0.10084 0.31678
{'drug_on_Drug on:attention_Attend RF'} 0.093561 0.10547 0.88713 120 0.37678 -0.11525 0.30238
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.0286 0.73701 1.4355
Group: Error
Name Estimate Lower Upper
{'Res Std'} 1.1744 1.0172 1.3559
Linear mixed-effects model fit by ML
Model information:
Number of observations 118
Fixed effects coefficients 4
Random effects coefficients 31
Covariance parameters 2
Formula:
gain_log ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
238.62 255.24 -113.31 226.62
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.8048 0.14524 -12.426 114 6.1881e-23 -2.0925 -1.5171
{'drug_on_Drug on' } 0.19251 0.042312 4.5498 114 1.3492e-05 0.10869 0.27633
{'attention_Attend RF' } -0.088546 0.042312 -2.0927 114 0.038593 -0.17236 -0.0047272
{'drug_on_Drug on:attention_Attend RF'} 0.023579 0.042629 0.55313 114 0.58126 -0.060868 0.10803
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 0.77318 0.58935 1.0144
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.4571 0.39407 0.53022
AUROC: Dopamine, delta-AUROC -0.02 +- 0.01, p = 0.094, Cohens's d = -0.28
Linear mixed-effects model fit by ML
Model information:
Number of observations 56
Fixed effects coefficients 4
Random effects coefficients 14
Covariance parameters 2
Formula:
rate ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
355.56 367.71 -171.78 343.56
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 15.922 2.6476 6.014 52 1.8299e-07 10.61 21.235
{'drug_on_Drug on' } -1.3651 0.44484 -3.0688 52 0.0034106 -2.2578 -0.47248
{'attention_Attend RF' } 2.5624 0.44484 5.7602 52 4.5931e-07 1.6697 3.455
{'drug_on_Drug on:attention_Attend RF'} -0.16397 0.44484 -0.36861 52 0.71391 -1.0566 0.72866
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 9.7654 6.6701 14.297
Group: Error
Name Estimate Lower Upper
{'Res Std'} 3.3289 2.688 4.1226
Linear mixed-effects model fit by ML
Model information:
Number of observations 56
Fixed effects coefficients 4
Random effects coefficients 14
Covariance parameters 2
Formula:
FF ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
232.04 244.19 -110.02 220.04
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 3.8125 0.39513 9.6487 52 3.5217e-13 3.0196 4.6054
{'drug_on_Drug on' } 0.2818 0.19273 1.4622 52 0.14971 -0.10494 0.66855
{'attention_Attend RF' } 0.30423 0.19273 1.5785 52 0.12051 -0.082515 0.69097
{'drug_on_Drug on:attention_Attend RF'} 0.21862 0.19273 1.1343 52 0.26185 -0.16812 0.60536
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.2907 0.79022 2.108
Group: Error
Name Estimate Lower Upper
{'Res Std'} 1.4423 1.1646 1.7862
Linear mixed-effects model fit by ML
Model information:
Number of observations 53
Fixed effects coefficients 4
Random effects coefficients 14
Covariance parameters 2
Formula:
gain_log ~ 1 + drug_on*attention + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
131.34 143.16 -59.668 119.34
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.5478 0.28729 -5.3876 49 2.0228e-06 -2.1251 -0.97044
{'drug_on_Drug on' } 0.24421 0.071712 3.4054 49 0.0013266 0.1001 0.38832
{'attention_Attend RF' } -0.14014 0.071712 -1.9542 49 0.056402 -0.28425 0.003974
{'drug_on_Drug on:attention_Attend RF'} 0.15056 0.07154 2.1046 49 0.040483 0.0067948 0.29433
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.0411 0.70072 1.5467
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.51583 0.41309 0.64411
AUROC: SCH23390, delta-AUROC -0.02 +- 0.03, p = 0.583, Cohens's d = -0.20
ncol = 2 * (length(label_drug));
[fH, fSet] = plotj_initFig('figNum', 2, 'width', 20, 'height', 12, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap .* [1 0.8];
[P, h_text] = deal( NaN(length(label_drug), nrow, 2, 2) ); % drugtype, actvity-type, attention, drug-offon,
data2plot = att_table.(datatype);
for idrug = 1:length(label_drug)
subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabelDisplacement', [0.05, 0.02]);
set(gca,'units','normalized')
text(100, 400, ['\bf ' label_drug{idrug}], 'FontSize', fSet.Fontsize_title, 'Interpreter', 'tex')
title(sprintf('%s', label_attention{iatt}), 'FontSize', fSet.Fontsize_title)
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_att = att_table.attention == label_attention{iatt};
tmp_data = [data2plot(unit2plot & att_table.drug_on==label_drug_onoff{1}) ...
data2plot(unit2plot & att_table.drug_on==label_drug_onoff{2})];
plotj_scatter(tmp_data, ...
'MarkerSize', markersize, ...
'markerFaceColor', fSet.colors(1,:), 'markerFaceAlpha', 0.5, ...
P(idrug,irow,iatt) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');
d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');
delta_data = tmp_data(:,2)-tmp_data(:,1);
p_string = get_significance_strings(P(idrug,irow,iatt), 'rounded', 0);
fprintf('%s: %s, attend %s: delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
label_attention{iatt}, ...
nanstd(delta_data)/sqrt(length(tmp_data)), ...
x_pos = get_value_range(tmp_x, x_text);
y_pos = get_value_range(tmp_y, y_text);
h_text(idrug,irow,iatt) = text(x_pos, y_pos, ...
sprintf('%s', p_string), ...
'Color', fSet.colors(1,:));
x_pos = get_value_range(tmp_x, 0.6);
y_pos = get_value_range(tmp_y, 0.2);
text(x_pos, y_pos, sprintf('(n=%d)', length(tmp_data)), 'Color', fSet.colors(1,:));
set(gca,'YScale', scale, 'XScale', scale)
% % set(gca,'', [1 10 100])
set(gca,'Xtick', [1 10 100], 'XTickLabel', [1 10 100])
set(gca,'Ytick', [1 10 100], 'YTickLabel', [1 10 100])
xlabel('Firing rate no drug (Hz)','FontSize', fSet.Fontsize_text)
ylabel('Firing rate drug (Hz)','FontSize', fSet.Fontsize_text)
xlabel('Fano Factor no drug','FontSize', fSet.Fontsize_text)
ylabel('Fano Factor drug','FontSize', fSet.Fontsize_text)
xlabel('Gain variability no drug','FontSize', fSet.Fontsize_text)
ylabel('Gain variability drug','FontSize', fSet.Fontsize_text)
end
rate: Dopamine, attend Attend RF: delta-rate 5.30 +- 0.87, p < 0.001, Cohens's d = 1.09
rate: Dopamine, attend Attend away: delta-rate 4.15 +- 0.84, p < 0.001, Cohens's d = 0.88
rate: SCH23390, attend Attend RF: delta-rate 3.06 +- 0.80, p = 0.002, Cohens's d = 1.03
rate: SCH23390, attend Attend away: delta-rate 2.40 +- 0.91, p = 0.017, Cohens's d = 0.71
FF: Dopamine, attend Attend RF: delta-FF -0.14 +- 0.36, p = 0.779, Cohens's d = -0.07
FF: Dopamine, attend Attend away: delta-FF 0.23 +- 0.21, p = 0.504, Cohens's d = 0.20
FF: SCH23390, attend Attend RF: delta-FF -1.00 +- 0.79, p = 0.296, Cohens's d = -0.34
FF: SCH23390, attend Attend away: delta-FF -0.13 +- 0.33, p = 0.903, Cohens's d = -0.10
gain_log: Dopamine, attend Attend RF: delta-gain_log -0.46 +- 0.14, p = 0.006, Cohens's d = -0.57
gain_log: Dopamine, attend Attend away: delta-gain_log -0.34 +- 0.11, p = 0.006, Cohens's d = -0.59
gain_log: SCH23390, attend Attend RF: delta-gain_log -0.82 +- 0.31, p = 0.009, Cohens's d = -0.82
gain_log: SCH23390, attend Attend away: delta-gain_log -0.18 +- 0.12, p = 0.244, Cohens's d = -0.15
[p_fdr, p_masked] = FDR(P, 0.05);
plotj_text_emphasise(h_text, p_masked, 'italic');
plotj_text_emphasise(h_text, p_masked, 'bold');
colors = get_colors('spikewidth');
binedges = 50:binsize:600;
assert(any(ismember(binedges,waveform_cutoff)), 'threshold needs to be included in binedges')
[fH, fSet] = plotj_initFig('width', 20, 'height', 12, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap .* [1 0.8];
subtightplot(nrow, ncol, 1, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
plotj_initAx(fSet, 'axlabel', 1, 'axlabelDisplacement', [0.07, 0.01]);
x = waveform_time-(7*6*5.4);
plot(x, unitlist.waveform(unitlist.unit_class=='Narrow',:)*-1, 'Color', [colors(1,:) 0.5])
plot(x, unitlist.waveform(unitlist.unit_class=='Broad',:)*-1, 'Color', [colors(2,:) 0.5])
xlabel(['Time from peak of spike (' plotj_symbol('mu') 's)'], 'FontSize', fSet.Fontsize_text)
ylabel(['Normalized voltage'], 'FontSize', fSet.Fontsize_text)
xlim(minmax(waveform_time)+[0 -200]-7*6*5.4)
subtightplot(nrow, ncol, 2, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
plotj_initAx(fSet, 'axlabel', 2, 'axlabelDisplacement', [0.07, 0.01]);
h = plotj_hist({unitlist.peak_to_trough_time(unitlist.unit_class=='Narrow'), unitlist.peak_to_trough_time(unitlist.unit_class=='Broad')}, ...
'histstyle', 'stairs',...
'LineWidth', fSet.LineWidth);
xlabel(['Peak to through time (' char(181) 's)'], 'FontSize', fSet.Fontsize_text)
ylabel('Number of cells', 'FontSize', fSet.Fontsize_text)
hleg = legend(h, {'Narrow','Broad'}, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'AutoUpdate', 'Off', 'Location', 'East');
hleg.Position = hleg.Position + [0.08 0 0 0];
% do hartigan's dip test to see whether there is a significant dip in this
% [dip, p_value, xlow,xup]=HartigansDipSignifTest(allSpikeWaveform.width(idx_unit),5000);
[dip, p_value, xlow, xup, boot_dip]=CalibratedHartigansDipSignifTest(unitlist.peak_to_trough_time, 10000);
p_string = get_significance_strings(p_value, 'round', 0);
h_text(1) = text(400, 13, {'Calibrated Hartigan''s dip test:', p_string}, 'FontSize', fSet.Fontsize_text);
% Note these excellent explanations of effect coded predictors:
% - https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-interpret-the-coefficients-of-an-effect-coded-variable-involved-in-an-interaction-in-a-regression-model/#effect-coded-predictors-interacted-with-a-continuous-covariate
% - https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faqwhat-is-effect-coding/
data2plot = {'rate', 'FF', 'gain_log'};
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
tmp_table = att_table(idx_unit,:);
for idp = 1:length(data2plot)
fprintf('\nDRUG: %s, DATATYPE: %s, UNIT SELECTION: %s, n=%d \n', label_drug{idrug}, data2plot{idp}, selectivity_criterium, length(unique(att_table.unit(idx_unit))))
% model fit across attention, drug and unit_class
mod_population.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention*drug_on*unit_class + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
mod_population.anova = anova(mod_population.lme);
% % hierarchical model fits
% mod.m1_intercept.lme = fitlme(tmp_table,[data2plot{idp} ' ~ 1 + (1|unit)'],'DummyVarCoding','effects');
% mod.m2_att.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + (1|unit)'],'DummyVarCoding','effects');
% mod.m3_drug.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + (1|unit)'],'DummyVarCoding','effects');
% mod.m4_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + (1|unit)'],'DummyVarCoding','effects');
% mod.m5_att_drug.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention*drug_on + (1|unit)'],'DummyVarCoding','effects');
% mod.m6_att_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + (1|unit)'],'DummyVarCoding','effects');
% mod.m7_drug_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + drug_on:unit_class + (1|unit)'],'DummyVarCoding','effects');
% mod.m8_att_drug_unitc.lme = fitlme(tmp_table,[data2plot{idp} ' ~ attention + drug_on + unit_class + attention:drug_on + attention:unit_class + drug_on:unit_class + attention:drug_on:unit_class + (1|unit)'],'DummyVarCoding','effects');
% % likelihood ratio tests of hierarchical model fits
% mod.m2_att.compare = compare(mod.m1_intercept.lme, mod.m2_att.lme, 'CheckNesting', true);
% mod.m3_drug.compare = compare(mod.m2_att.lme, mod.m3_drug.lme, 'CheckNesting', true);
% mod.m4_unitc.compare = compare(mod.m3_drug.lme, mod.m4_unitc.lme, 'CheckNesting', true);
% mod.m5_att_drug.compare = compare(mod.m4_unitc.lme, mod.m5_att_drug.lme, 'CheckNesting', true);
% mod.m6_att_unitc.compare = compare(mod.m5_att_drug.lme, mod.m6_att_unitc.lme, 'CheckNesting', true);
% mod.m7_drug_unitc.compare = compare(mod.m6_att_unitc.lme, mod.m7_drug_unitc.lme, 'CheckNesting', true);
% mod.m8_att_drug_unitc.compare = compare(mod.m7_drug_unitc.lme, mod.m8_att_drug_unitc.lme, 'CheckNesting', true);
% print stats for each group
disp(grpstats(tmp_table,{'attention','drug_on','unit_class'},{'min','max','mean','median'},'DataVars',data2plot{idp}))
end
DRUG: Dopamine, DATATYPE: rate, UNIT SELECTION: att&dru, n=31
attention drug_on unit_class GroupCount min_rate max_rate mean_rate median_rate
___________ ________ __________ __________ ________ ________ _________ ___________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 14 0.79125 47.567 16.942 15.312
Attend RF_Drug on_Broad Attend RF Drug on Broad 17 3.4929 69.675 21.011 10.703
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 14 1.5655 52.707 20.546 21.285
Attend RF_Drug off_Broad Attend RF Drug off Broad 17 4.5532 80.411 27.717 15.746
Attend away_Drug on_Narrow Attend away Drug on Narrow 14 1.9689 34.203 12.246 9.0068
Attend away_Drug on_Broad Attend away Drug on Broad 17 4.1932 54.01 16.105 8.9194
Attend away_Drug off_Narrow Attend away Drug off Narrow 14 4.2369 38.811 15.463 13.631
Attend away_Drug off_Broad Attend away Drug off Broad 17 5.637 60.66 21.026 11.448
Linear mixed-effects model fit by ML
Model information:
Number of observations 124
Fixed effects coefficients 8
Random effects coefficients 31
Covariance parameters 2
Formula:
rate ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
852.27 880.47 -416.13 832.27
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 18.882 2.724 6.9317 116 2.4955e-10 13.487 24.277
{'drug_on_Drug on' } -2.3057 0.38341 -6.0137 116 2.1581e-08 -3.0651 -1.5463
{'attention_Attend RF' } 2.6719 0.38341 6.9688 116 2.0728e-10 1.9125 3.4313
{'unit_class_Narrow' } -2.5828 2.724 -0.94818 116 0.34501 -7.9781 2.8124
{'drug_on_Drug on:attention_Attend RF' } -0.2715 0.38341 -0.7081 116 0.4803 -1.0309 0.4879
{'drug_on_Drug on:unit_class_Narrow' } 0.60072 0.38341 1.5668 116 0.11989 -0.15868 1.3601
{'attention_Attend RF:unit_class_Narrow' } -0.22728 0.38341 -0.59279 116 0.55448 -0.98668 0.53212
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} 0.1747 0.38341 0.45564 116 0.6495 -0.5847 0.9341
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 14.945 11.593 19.266
Group: Error
Name Estimate Lower Upper
{'Res Std'} 4.2495 3.6806 4.9062
DRUG: Dopamine, DATATYPE: FF, UNIT SELECTION: att&dru, n=31
attention drug_on unit_class GroupCount min_FF max_FF mean_FF median_FF
___________ ________ __________ __________ ______ ______ _______ _________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 14 2.0911 7.1985 3.8655 3.0195
Attend RF_Drug on_Broad Attend RF Drug on Broad 17 1.2696 6.3605 3.5064 3.3962
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 14 1.4285 5.2381 2.8399 2.3452
Attend RF_Drug off_Broad Attend RF Drug off Broad 17 2.0925 10.152 4.0938 3.2168
Attend away_Drug on_Narrow Attend away Drug on Narrow 14 1.7752 5.4315 3.213 2.8758
Attend away_Drug on_Broad Attend away Drug on Broad 17 2.0443 5.7579 3.3087 2.6984
Attend away_Drug off_Narrow Attend away Drug off Narrow 14 1.8161 6.673 3.5405 3.2446
Attend away_Drug off_Broad Attend away Drug off Broad 17 1.4416 8.299 3.4643 2.8067
Linear mixed-effects model fit by ML
Model information:
Number of observations 124
Fixed effects coefficients 8
Random effects coefficients 31
Covariance parameters 2
Formula:
FF ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
446.19 474.39 -213.09 426.19
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 3.479 0.21274 16.353 116 6.766e-32 3.0576 3.9004
{'drug_on_Drug on' } -0.0056168 0.10107 -0.055573 116 0.95578 -0.2058 0.19457
{'attention_Attend RF' } 0.097384 0.10107 0.96352 116 0.33729 -0.1028 0.29757
{'unit_class_Narrow' } -0.11428 0.21274 -0.53717 116 0.59218 -0.53564 0.30708
{'drug_on_Drug on:attention_Attend RF' } 0.11515 0.10107 1.1393 116 0.25692 -0.085032 0.31533
{'drug_on_Drug on:unit_class_Narrow' } 0.18014 0.10107 1.7823 116 0.077315 -0.020044 0.38032
{'attention_Attend RF:unit_class_Narrow' } -0.1094 0.10107 -1.0824 116 0.28132 -0.30958 0.090784
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} 0.22309 0.10107 2.2073 116 0.029263 0.022908 0.42327
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.0374 0.75014 1.4346
Group: Error
Name Estimate Lower Upper
{'Res Std'} 1.1202 0.97024 1.2933
DRUG: Dopamine, DATATYPE: gain_log, UNIT SELECTION: att&dru, n=31
attention drug_on unit_class GroupCount min_gain_log max_gain_log mean_gain_log median_gain_log
___________ ________ __________ __________ ____________ ____________ _____________ _______________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 14 -3.1554 1.0697 -1.4737 -1.4938
Attend RF_Drug on_Broad Attend RF Drug on Broad 17 -4.2117 -0.46886 -1.8043 -1.9479
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 14 -3.2009 -0.0074051 -2.115 -2.3693
Attend RF_Drug off_Broad Attend RF Drug off Broad 17 -3.6859 -0.82216 -2.0801 -2.031
Attend away_Drug on_Narrow Attend away Drug on Narrow 14 -2.9678 0.86773 -1.2674 -1.4614
Attend away_Drug on_Broad Attend away Drug on Broad 17 -3.4723 -0.53745 -1.6883 -1.6183
Attend away_Drug off_Narrow Attend away Drug off Narrow 14 -3.2715 -0.11799 -1.6461 -1.8262
Attend away_Drug off_Broad Attend away Drug off Broad 17 -3.5353 -0.76188 -2.0821 -2.0195
Linear mixed-effects model fit by ML
Model information:
Number of observations 118
Fixed effects coefficients 8
Random effects coefficients 31
Covariance parameters 2
Formula:
gain_log ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
240.21 267.92 -110.1 220.21
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.7965 0.14453 -12.43 110 1.0898e-22 -2.0829 -1.5101
{'drug_on_Drug on' } 0.20057 0.0414 4.8446 110 4.1833e-06 0.11852 0.28261
{'attention_Attend RF' } -0.099717 0.0414 -2.4086 110 0.017675 -0.18176 -0.017672
{'unit_class_Narrow' } 0.11883 0.14453 0.82221 110 0.41274 -0.16759 0.40526
{'drug_on_Drug on:attention_Attend RF' } 0.033293 0.041782 0.79683 110 0.42726 -0.049509 0.11609
{'drug_on_Drug on:unit_class_Narrow' } 0.034814 0.0414 0.8409 110 0.40223 -0.047232 0.11686
{'attention_Attend RF:unit_class_Narrow' } -0.06958 0.0414 -1.6807 110 0.095665 -0.15163 0.012465
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} 0.06441 0.041782 1.5416 110 0.12605 -0.018392 0.14721
Random effects covariance parameters (95% CIs):
Group: unit (31 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 0.76695 0.58502 1.0055
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.44218 0.38117 0.51295
DRUG: SCH23390, DATATYPE: rate, UNIT SELECTION: att&dru, n=14
attention drug_on unit_class GroupCount min_rate max_rate mean_rate median_rate
___________ ________ __________ __________ ________ ________ _________ ___________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 3 10.527 38.003 19.858 11.043
Attend RF_Drug on_Broad Attend RF Drug on Broad 11 1.8333 31.818 16.164 16.361
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 3 12.542 43.356 23.049 13.248
Attend RF_Drug off_Broad Attend RF Drug off Broad 11 2.7802 37.4 19.186 17.697
Attend away_Drug on_Narrow Attend away Drug on Narrow 3 4.9403 21.982 11.365 7.1728
Attend away_Drug on_Broad Attend away Drug on Broad 11 1.0897 26.58 12.375 8.2434
Attend away_Drug off_Narrow Attend away Drug off Narrow 3 6.5644 24.866 12.802 6.976
Attend away_Drug off_Broad Attend away Drug off Broad 11 1.6288 34.826 15.041 11.932
Linear mixed-effects model fit by ML
Model information:
Number of observations 56
Fixed effects coefficients 8
Random effects coefficients 14
Covariance parameters 2
Formula:
rate ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
356.62 376.87 -168.31 336.62
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 16.23 3.2229 5.0358 48 7.1363e-06 9.7498 22.71
{'drug_on_Drug on' } -1.2895 0.49925 -2.5828 48 0.012903 -2.2933 -0.28567
{'attention_Attend RF' } 3.3342 0.49925 6.6784 48 2.2901e-08 2.3304 4.338
{'unit_class_Narrow' } 0.5384 3.2229 0.16705 48 0.86803 -5.9418 7.0186
{'drug_on_Drug on:attention_Attend RF' } -0.26378 0.49925 -0.52836 48 0.59968 -1.2676 0.74002
{'drug_on_Drug on:unit_class_Narrow' } 0.13238 0.49925 0.26515 48 0.79203 -0.87143 1.1362
{'attention_Attend RF:unit_class_Narrow' } 1.3506 0.49925 2.7054 48 0.0094162 0.34684 2.3545
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} -0.17466 0.49925 -0.34986 48 0.72798 -1.1785 0.82914
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 9.7769 6.6891 14.29
Group: Error
Name Estimate Lower Upper
{'Res Std'} 3.066 2.4757 3.797
DRUG: SCH23390, DATATYPE: FF, UNIT SELECTION: att&dru, n=14
attention drug_on unit_class GroupCount min_FF max_FF mean_FF median_FF
___________ ________ __________ __________ ______ ______ _______ _________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 3 4.8438 6.9749 5.6723 5.1983
Attend RF_Drug on_Broad Attend RF Drug on Broad 11 1.2078 14.527 4.3294 3.6886
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 3 4.7386 7.6174 5.9872 5.6056
Attend RF_Drug off_Broad Attend RF Drug off Broad 11 1.2523 4.7912 2.9697 2.8743
Attend away_Drug on_Narrow Attend away Drug on Narrow 3 3.2601 4.786 4.1195 4.3125
Attend away_Drug on_Broad Attend away Drug on Broad 11 2.1601 7.0603 3.422 2.7141
Attend away_Drug off_Narrow Attend away Drug off Narrow 3 4.0334 5.6038 4.7933 4.7427
Attend away_Drug off_Broad Attend away Drug off Broad 11 1.8187 4.2039 3.0774 3.1568
Linear mixed-effects model fit by ML
Model information:
Number of observations 56
Fixed effects coefficients 8
Random effects coefficients 14
Covariance parameters 2
Formula:
FF ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
233.16 253.42 -106.58 213.16
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } 4.2964 0.42499 10.109 48 1.7726e-13 3.4419 5.1509
{'drug_on_Drug on' } 0.089451 0.22558 0.39654 48 0.69347 -0.36411 0.54301
{'attention_Attend RF' } 0.4433 0.22558 1.9652 48 0.055197 -0.010258 0.89686
{'unit_class_Narrow' } 0.84672 0.42499 1.9923 48 0.05204 -0.0077845 1.7012
{'drug_on_Drug on:attention_Attend RF' } 0.17175 0.22558 0.76135 48 0.45017 -0.28181 0.62531
{'drug_on_Drug on:unit_class_Narrow' } -0.33662 0.22558 -1.4922 48 0.14218 -0.79018 0.11694
{'attention_Attend RF:unit_class_Narrow' } 0.24338 0.22558 1.0789 48 0.28602 -0.21018 0.69694
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} -0.082033 0.22558 -0.36365 48 0.71771 -0.53559 0.37153
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.106 0.65592 1.8649
Group: Error
Name Estimate Lower Upper
{'Res Std'} 1.3853 1.1186 1.7156
DRUG: SCH23390, DATATYPE: gain_log, UNIT SELECTION: att&dru, n=14
attention drug_on unit_class GroupCount min_gain_log max_gain_log mean_gain_log median_gain_log
___________ ________ __________ __________ ____________ ____________ _____________ _______________
Attend RF_Drug on_Narrow Attend RF Drug on Narrow 3 -2.0016 -0.65222 -1.2007 -0.94844
Attend RF_Drug on_Broad Attend RF Drug on Broad 11 -3.5623 1.1763 -1.121 -1.6503
Attend RF_Drug off_Narrow Attend RF Drug off Narrow 3 -1.9344 -0.89988 -1.3305 -1.1571
Attend RF_Drug off_Broad Attend RF Drug off Broad 11 -4.0364 -0.22889 -2.2878 -2.6556
Attend away_Drug on_Narrow Attend away Drug on Narrow 3 -1.9583 -0.54467 -1.0418 -0.62238
Attend away_Drug on_Broad Attend away Drug on Broad 11 -3.1992 1.2952 -1.3882 -1.477
Attend away_Drug off_Narrow Attend away Drug off Narrow 3 -1.679 -0.38222 -0.92504 -0.71391
Attend away_Drug off_Broad Attend away Drug off Broad 11 -3.2254 0.82684 -1.5237 -1.3357
Linear mixed-effects model fit by ML
Model information:
Number of observations 53
Fixed effects coefficients 8
Random effects coefficients 14
Covariance parameters 2
Formula:
gain_log ~ 1 + drug_on*attention + drug_on*unit_class + attention*unit_class + drug_on:attention:unit_class + (1 | unit)
Model fit statistics:
AIC BIC LogLikelihood Deviance
134.64 154.34 -57.318 114.64
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.3918 0.3414 -4.0767 45 0.00018358 -2.0794 -0.70419
{'drug_on_Drug on' } 0.16056 0.080691 1.9898 45 0.052716 -0.0019651 0.32308
{'attention_Attend RF' } -0.13652 0.080691 -1.6918 45 0.09759 -0.29904 0.026004
{'unit_class_Narrow' } 0.2673 0.3414 0.78294 45 0.43776 -0.42033 0.95492
{'drug_on_Drug on:attention_Attend RF' } 0.12026 0.080615 1.4918 45 0.14272 -0.042104 0.28263
{'drug_on_Drug on:unit_class_Narrow' } -0.15731 0.080691 -1.9496 45 0.057475 -0.31983 0.0052072
{'attention_Attend RF:unit_class_Narrow' } -0.0045813 0.080691 -0.056775 45 0.95498 -0.1671 0.15794
{'drug_on_Drug on:attention_Attend RF:unit_class_Narrow'} -0.058643 0.080615 -0.72744 45 0.47073 -0.22101 0.10372
Random effects covariance parameters (95% CIs):
Group: unit (14 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 1.0187 0.68687 1.5108
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.48979 0.39227 0.61156
data2plot = {'rate','FF', 'gain_log'};
[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');
colors_violin = reshape(colors_violin, [4,3]);
labels_violin = reshape(labels_violin, [4,1]);
labels_violin_short = reshape(labels_violin_short, [4,1]);
% plotting and post-hoc tests
ncol = num_drug*length(label_unitclass); % data2plot x unit
nrow = length(data2plot);
% use subplot indices to combine violin and post-hoc scatter plots
idx_axlabel = [1 NaN nrow+1 NaN 2 NaN nrow+2 NaN 3 NaN nrow+3 NaN];
[fH, fSet] = plotj_initFig('figNum', 3, 'width', 20, 'height', 15, 'Journal', plot_conventions);
fSet.subplotGap = [0.1, 0.06];
for idp = 1:length(data2plot)
axpos = zeros(length(label_unitclass),4);
fprintf('Plot data for drug: %s \n', label_drug{idrug})
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
fprintf('\t Found %d units \n', length(unique(att_table.unit(idx_unit))))
% plot, violinplot for each unitclass
for iunitc = 1:length(label_unitclass) % for each unit_class
subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabel', idx_axlabel(iplot), 'axlabelDisplacement', [0.05, 0.02]);
axH = plotj_initAx(fSet);
axH.Position = axH.Position - [0.08 0 0 0];
h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);
idx_unit_class = att_table.unit_class==label_unitclass(iunitc);
num_unit_plot = length(unique(att_table.unit(idx_unit & idx_unit_class)));
v_att_table = att_table(idx_unit & idx_unit_class, ismember(att_table.Properties.VariableNames, {data2plot{idp},'drug_on','attention','unit','attention_drug_on'}));
fprintf('\t Found %d independent units in unit_class: %s \n', length(unique(v_att_table.unit)), label_unitclass{iunitc})
v = violinplot(v_att_table.(data2plot{idp}), v_att_table.attention_drug_on, 'ShowMean', true);
XTickLabel = get(gca,'XTickLabel');
v(iv).ScatterPlot.SizeData = 10;
idx = strcmpi(XTickLabel(iv), labels_violin);
v(iv).ViolinColor = colors_violin(idx,:);
% plot lines for each unit
[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));
if ~any(isnan(v_att_table.(data2plot{idp})))
% works if data is complete
data_x(:,iv) = v(iv).ScatterPlot.XData;
data_y(:,iv) = v(iv).ScatterPlot.YData;
% catch if nan values are present
idx = strcmpi(XTickLabel(iv), labels_violin);
idx_trial = v_att_table.attention_drug_on==labels_violin{idx};
data_y(:,iv) = v_att_table.(data2plot{idp})(idx_trial);
plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);
lme_rate = fitlme(v_att_table,[data2plot{idp} ' ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
% fprintf('Model fit: %s, unit: %s\n', label_drug{idrug}, label_unittype{iunitc})
coefficients = lme_rate.Coefficients;
betas = coefficients.Estimate(2:end);
betas_se = coefficients.SE(2:end);
p_val = coefficients.pValue(2:end);
% p_string = get_significance_strings(p_val, 'rounded', 0);
p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});
p_string{end+1} = sprintf('n=%d',num_unit_plot);
% % full betastring: Estimate, SE, p-value
% sprintf('drug: %s = %1.2f%s%1.2f, %s', '\beta', betas(1), '\pm', betas_se(1), p_string{1}), ...
% sprintf('att: %s = %1.2f%s%1.2f, %s', '\beta', betas(2), '\pm', betas_se(2), p_string{2}), ...
% sprintf('att %s drug: %s = %1.2f%s%1.2f, %s', '\times', '\beta', betas(3), '\pm', betas_se(3), p_string{3}), ...
% sprintf('n = %d',length(unique(v_att_table.unit)))};
% % part betastring: Estimate, SE, p-value
% sprintf('drug: %s = %1.2f, %s', '\beta', betas(1), p_string{1}), ...
% sprintf('att: %s = %1.2f, %s', '\beta', betas(2), p_string{2}), ...
% sprintf('att %s drug: %s = %1.2f, %s', '\times', '\beta', betas(3), p_string{3}), ...
% sprintf('n = %d',length(unique(v_att_table.unit)))};
set(gca, 'ylim', [0 1000])
set(gca, 'ylim', YLIM .* [1 1.5])
x_pos = get_value_range(tmp_x, 1);
y_pos = get_value_range(tmp_y, y_text);
h_text = text(x_pos, y_pos, p_string, ...
'HorizontalAlignment','right',...
'FontSize',fSet.Fontsize_text_in);
ylabel('Firing rate (Hz)', 'FontSize', fSet.Fontsize_text)
ylabel('Fano Factor', 'FontSize', fSet.Fontsize_text)
ylabel('\sigma_{G}^{2}', 'FontSize', fSet.Fontsize_text)
ylabel('log(\sigma_{G}^{2})', 'FontSize', fSet.Fontsize_text)
if (idp<length(data2plot)) || (idrug>1) || (iunitc>1)
set(gca,'XTickLabelRotation',-25)
% [p_fdr, p_masked] = FDR(p_val, 0.05);
p_masked = [p_masked ; false];
plotj_text_emphasise(h_text, 1, 'italic', p_masked);
plotj_text_emphasise(h_text, 1, 'bold', p_masked);
end
Plot data for drug: Dopamine
Found 14 independent units in unit_class: Narrow
Found 17 independent units in unit_class: Broad
Plot data for drug: SCH23390
Found 3 independent units in unit_class: Narrow
Found 11 independent units in unit_class: Broad
Plot data for drug: Dopamine
Found 14 independent units in unit_class: Narrow
Found 17 independent units in unit_class: Broad
Plot data for drug: SCH23390
Found 3 independent units in unit_class: Narrow
Found 11 independent units in unit_class: Broad
Plot data for drug: Dopamine
Found 14 independent units in unit_class: Narrow
Found 17 independent units in unit_class: Broad
Plot data for drug: SCH23390
Found 3 independent units in unit_class: Narrow
Found 11 independent units in unit_class: Broad
ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized','clipping','off');
xpos = (1/(ncol)) * xpos_idx(idrug);
ht = text(xpos, 0.95, [label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title * 1.2, 'HorizontalAlignment', 'center');
ncol = 2 * (length(label_drug));
colors_scatter = get_colors('spikewidth');
idx_subplot = [1 2 ncol+1 ncol+2 2*ncol+1 2*ncol+2; 3 4 ncol+3 ncol+4 2*ncol+3 2*ncol+4] ;
idx_axlabel = [1 NaN 2 NaN 3 NaN; 4 NaN 5 NaN 6 NaN];
[fH, fSet] = plotj_initFig('figNum', 2, 'width', 20, 'height', 16, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap .* [1 0.8];
[P, h_text] = deal( NaN(length(label_drug), nrow, 2, 2) ); % drugtype, actvity-type, attention, drug-offon,
for idrug = 1:length(label_drug)
data2plot = att_table.(datatype);
subtightplot(nrow, ncol, idx_subplot(idrug, iplot), fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
axH = plotj_initAx(fSet, 'axlabel', idx_axlabel(idrug, iplot), 'axlabelDisplacement', [0.05, 0.02]);
set(gca,'units','normalized')
text(100, 400, ['\bf ' label_drug{idrug}], 'FontSize', fSet.Fontsize_title, 'Interpreter', 'tex')
title(sprintf('Attend %s', label_attention{iatt}), 'FontSize', fSet.Fontsize_title)
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_att = att_table.attention == label_attention{iatt};
idx_unit_class = att_table.unit_class == label_unitclass{iunitc};
tmp_data = [data2plot(unit2plot & att_table.drug_on==label_drug_onoff{1}) ...
data2plot(unit2plot & att_table.drug_on==label_drug_onoff{2})];
plotj_scatter(tmp_data, ...
'markerStyle', markerstyle(iunitc), 'MarkerSize', markersize, ...
'markerFaceColor', colors_scatter(iunitc,:), 'markerFaceAlpha', 0.5, ...
'markerEdgeColor', colors_scatter(iunitc,:), 'markerEdgeAlpha', 0.5, ...
P(idrug,irow,iatt,iunitc) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');
d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');
delta_data = tmp_data(:,2)-tmp_data(:,1);
p_string = get_significance_strings(P(idrug,irow,iatt,iunitc), 'rounded', 0);
fprintf('%s: %s, %s, attend %s: delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
label_unitclass{iunitc}, ...
label_attention{iatt}, ...
std(delta_data)/sqrt(length(tmp_data)), ...
x_pos = get_value_range(tmp_x, x_text);
y_pos = get_value_range(tmp_y, y_text(1)-y_text(2)*(iunitc-1));
h_text(idrug,irow,iatt,iunitc) = text(x_pos, y_pos, ...
sprintf('%s', p_string), ...
'Color', colors_scatter(iunitc,:));
x_pos = get_value_range(tmp_x, 0.6);
y_pos = get_value_range(tmp_y, 0.1+0.1*(iunitc-1));
text(x_pos, y_pos, sprintf('%s (n=%d)', label_unitclass{iunitc}, length(tmp_data)), 'Color', colors_scatter(iunitc,:));
set(gca,'YScale', scale, 'XScale', scale)
% % set(gca,'', [1 10 100])
if strcmpi(datatype, 'rate')
set(gca,'Xtick', [1 10 100], 'XTickLabel', [1 10 100])
set(gca,'Ytick', [1 10 100], 'YTickLabel', [1 10 100])
xlabel('Firing rate no drug (Hz)','FontSize', fSet.Fontsize_text)
ylabel('Firing rate drug (Hz)','FontSize', fSet.Fontsize_text)
elseif strcmp(datatype, 'FF')
xlabel('Fano Factor no drug','FontSize', fSet.Fontsize_text)
ylabel('Fano Factor drug','FontSize', fSet.Fontsize_text)
elseif strcmp(datatype, 'gain_log')
xlabel('Gain no drug','FontSize', fSet.Fontsize_text)
ylabel('Gain drug','FontSize', fSet.Fontsize_text)
end
rate: Dopamine, Narrow, attend Attend RF: delta-rate 3.60 +- 1.12, p = 0.009, Cohens's d = 0.86
rate: Dopamine, Broad, attend Attend RF: delta-rate 6.71 +- 1.22, p < 0.001, Cohens's d = 1.33
rate: Dopamine, Narrow, attend Attend away: delta-rate 3.22 +- 1.06, p = 0.009, Cohens's d = 0.81
rate: Dopamine, Broad, attend Attend away: delta-rate 4.92 +- 1.26, p < 0.001, Cohens's d = 0.95
FF: Dopamine, Narrow, attend Attend RF: delta-FF -1.03 +- 0.58, p = 0.119, Cohens's d = -0.47
FF: Dopamine, Broad, attend Attend RF: delta-FF 0.59 +- 0.39, p = 0.174, Cohens's d = 0.37
FF: Dopamine, Narrow, attend Attend away: delta-FF 0.33 +- 0.30, p = 0.326, Cohens's d = 0.30
FF: Dopamine, Broad, attend Attend away: delta-FF 0.16 +- 0.31, p = 0.927, Cohens's d = 0.12
gain_log: Dopamine, Narrow, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.032, Cohens's d = -0.64
gain_log: Dopamine, Broad, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.093, Cohens's d = -0.54
gain_log: Dopamine, Narrow, attend Attend away: delta-gain_log NaN +- NaN, p = 0.204, Cohens's d = -0.55
gain_log: Dopamine, Broad, attend Attend away: delta-gain_log -0.39 +- 0.15, p = 0.011, Cohens's d = -0.65
rate: SCH23390, Narrow, attend Attend RF: delta-rate 3.19 +- 1.14, p = 0.250, Cohens's d = 1.62
rate: SCH23390, Broad, attend Attend RF: delta-rate 3.02 +- 0.99, p = 0.010, Cohens's d = 0.92
rate: SCH23390, Narrow, attend Attend away: delta-rate 1.44 +- 0.89, p = 0.500, Cohens's d = 0.93
rate: SCH23390, Broad, attend Attend away: delta-rate 2.67 +- 1.13, p = 0.042, Cohens's d = 0.71
FF: SCH23390, Narrow, attend Attend RF: delta-FF 0.31 +- 0.39, p = 0.500, Cohens's d = 0.47
FF: SCH23390, Broad, attend Attend RF: delta-FF -1.36 +- 0.98, p = 0.206, Cohens's d = -0.42
FF: SCH23390, Narrow, attend Attend away: delta-FF 0.67 +- 0.39, p = 0.500, Cohens's d = 1.00
FF: SCH23390, Broad, attend Attend away: delta-FF -0.34 +- 0.39, p = 0.638, Cohens's d = -0.26
gain_log: SCH23390, Narrow, attend Attend RF: delta-gain_log -0.13 +- 0.10, p = 0.500, Cohens's d = -0.76
gain_log: SCH23390, Broad, attend Attend RF: delta-gain_log NaN +- NaN, p = 0.020, Cohens's d = -0.94
gain_log: SCH23390, Narrow, attend Attend away: delta-gain_log 0.12 +- 0.11, p = 0.500, Cohens's d = 0.62
gain_log: SCH23390, Broad, attend Attend away: delta-gain_log NaN +- NaN, p = 0.160, Cohens's d = -0.27
[p_fdr, p_masked] = FDR(P, 0.05);
plotj_text_emphasise(h_text, p_masked, 'italic');
plotj_text_emphasise(h_text, p_masked, 'bold');
[colors, labels] = get_colors('att_drug');
[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');
colors_violin = reshape(colors_violin, [4,3]);
labels_violin = reshape(labels_violin, [4,1]);
ncol = length(label_drug);
[fH, fSet] = plotj_initFig('width',20,'height',15,'Journal',plot_conventions);
for idrug = 1:length(label_drug)
for iunitc = 1:length(unique(unitlist.unit_class))
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);
h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', idrug);
h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);
gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
variable_names = gv.Properties.VariableNames;
for iatt_drug = 1:length(labels_violin)
mus = [gv.(labels_violin{iatt_drug}).mus];
s2s = [gv.(labels_violin{iatt_drug}).s2s];
rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);
gv_pop = GainVariability;
m2v = gv_pop.pred_variance(xval_gain);
h(iatt) = scatter(mus,s2s,fSet.MarkerSize*3,colors_violin(iatt_drug,:),'filled','MarkerFaceAlpha',0.5);
plot(m2v.mean,m2v.variance,'Color',colors_violin(iatt_drug,:),'linew',fSet.LineWidth);
xlim([0.8 max(xval_gain)])
ylim([0.8 max(xval_gain)])
plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)
set(gca,'xscale','log','yscale','log')
% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)
if idrug==length(label_drug)
xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)
ylabel({[label_drug{idrug} label_drug_ext{idrug}],'Variance (spikes^{2})'},'FontSize',fSet.Fontsize_text)
pos = get(h_ax, 'Position');
h_ax_inset = axes('Position', [pos(1)+0.22 pos(2)+0.03 0.16 0.15]) ; % inset
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);
num_unit_plot = length(unique(att_table.unit(idx_unit)));
fprintf('\t Found %d units \n', num_unit_plot)
v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {'gain_log','drug_on','attention','unit','attention_drug_on'}));
v = violinplot(v_att_table.gain_log, v_att_table.attention_drug_on, 'ShowMean', true);
XTickLabel = get(gca,'XTickLabel');
v(iv).ScatterPlot.SizeData = 10;
idx = strcmpi(XTickLabel(iv), labels_violin);
v(iv).ViolinColor = colors_violin(idx,:);
h_ax_inset.XAxis.Visible = 'off';
% plot lines for each unit
[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));
if ~any(isnan(v_att_table.gain_log))
% works if data is complete
data_x(:,iv) = v(iv).ScatterPlot.XData;
data_y(:,iv) = v(iv).ScatterPlot.YData;
% catch if nan values are present
idx = strcmpi(XTickLabel(iv), labels_violin);
idx_trial = v_att_table.attention_drug_on==labels_violin{idx};
data_y(:,iv) = v_att_table.gain_log(idx_trial);
plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);
lme_population = fitlme(v_att_table,['gain_log ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
disp(lme_population.Coefficients)
coefficients = lme_population.Coefficients;
betas = coefficients.Estimate(2:end);
betas_se = coefficients.SE(2:end);
p_val = coefficients.pValue(2:end);
% p_string = get_significance_strings(p_val, 'rounded', 0);
p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});
p_string{end+1} = ['n = ' num2str(num_unit_plot)];
set(gca, 'ylim', YLIM .* [1 1.25])
x_pos = get_value_range(tmp_x, 0.95);
y_pos = get_value_range(tmp_y, 1.3);
h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right');
ylabel('log(\sigma^{2}_G )', 'FontSize', fSet.Fontsize_text)
% [p_fdr, p_masked] = FDR(p_val, 0.05);
p_masked = [p_masked ; false];
plotj_text_emphasise(h_text, 1, 'italic', p_masked);
plotj_text_emphasise(h_text, 1, 'bold', p_masked);
end
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.6769 0.24109 -6.9553 47 9.5535e-09 -2.1619 -1.1919
{'drug_on_Drug on' } 0.23563 0.074521 3.1619 47 0.002744 0.085713 0.38555
{'attention_Attend RF' } -0.16931 0.074521 -2.272 47 0.027706 -0.31923 -0.019394
{'drug_on_Drug on:attention_Attend RF'} 0.097133 0.075694 1.2832 47 0.20571 -0.055145 0.24941
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.9154 0.1718 -11.149 63 1.5163e-16 -2.2587 -1.572
{'drug_on_Drug on' } 0.16573 0.045026 3.6809 63 0.00048392 0.075757 0.25571
{'attention_Attend RF' } -0.030159 0.045026 -0.66981 63 0.50543 -0.12014 0.059818
{'drug_on_Drug on:attention_Attend RF'} -0.031138 0.045026 -0.69156 63 0.49176 -0.12111 0.058839
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.1245 0.31895 -3.5257 8 0.0077812 -1.86 -0.38901
{'drug_on_Drug on' } 0.0032421 0.029576 0.10962 8 0.91541 -0.064961 0.071446
{'attention_Attend RF' } -0.1411 0.029576 -4.7706 8 0.0014075 -0.2093 -0.072895
{'drug_on_Drug on:attention_Attend RF'} 0.06162 0.029576 2.0834 8 0.070737 -0.0065831 0.12982
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.6587 0.34668 -4.7845 37 2.7399e-05 -2.3611 -0.95625
{'drug_on_Drug on' } 0.31791 0.088241 3.6027 37 0.00092093 0.13911 0.4967
{'attention_Attend RF' } -0.1319 0.088241 -1.4948 37 0.14346 -0.31069 0.046895
{'drug_on_Drug on:attention_Attend RF'} 0.17933 0.087884 2.0405 37 0.048482 0.0012596 0.3574
% datatype for dose-response curve
ncol = length(label_drug) ;
[fH, fSet] = plotj_initFig('width', 15, 'height', 8, 'Journal',plot_conventions);
v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);
ylabel2use = 'Drug modulation index';
ylim2use = {[-0.5 0.5],[-0.4 0.2]};
idx_cond = att_table.attention=='Attend RF';
tmp_data = zeros(length(find(idx_cond))/2,1);
tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');
tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');
data2plot = diff(tmp_data,1,2);
ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];
ylim2use = {[-0.25 0.3],[-0.25 0.3]};
v_att_table = unstack(att_table, 'gain_log', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
% difference in attentional modulation of gain variabiliy
data2plot1 = (v_att_table.("Attend RF, Drug on") - v_att_table.("Attend away, Drug on")) ./ (v_att_table.("Attend RF, Drug on") + v_att_table.("Attend away, Drug on"));
data2plot2 = (v_att_table.("Attend RF, Drug off") - v_att_table.("Attend away, Drug off")) ./ (v_att_table.("Attend RF, Drug off") + v_att_table.("Attend away, Drug off"));
data2plot = data2plot1 - data2plot2;
for idrug = 1:length(label_drug)
idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});
ax(idrug) = subtightplot(nrow, ncol, idrug, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', idrug);
title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
xlabel_event = 'Ejection Current';
tmp_data = [unitlist.EjectCurrent(idx_unit) data2plot(idx_unit)];
hscat = plotj_scatter(tmp_data, ...
'markerStyle', markerstyle, 'markerSize', markersize, ...
'markerFaceColor', 'k', 'markerFaceAlpha', 0.5, ...
text_legend{1} = ['(n=' num2str(length(find(idx_unit))) ')'];
data_table = table(unitlist.EjectCurrent_centered(idx_unit), data2plot(idx_unit), 'VariableNames', {'x', 'y'});
[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);
predict_data.x = predict_data.x + mean(unique(unitlist.EjectCurrent));
filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criterium));
filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criterium));
writetable(data_table, filename)
% fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)
if any(strcmp('fit_mean', predict_data.Properties.VariableNames))
h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);
cond1 = exist([filename(1:end-4) '_R.csv'], 'file');
cond2 = strcmpi(model,'quadratic');
utable = readtable([filename(1:end-4) '_R.csv']);
pstring_chi = [pstring_chi ' (U^{+})'];
pstring_chi = [pstring_chi ' (U^{-})'];
betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];
% model_color = colors(2,:);
betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];
xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)
ylabel(ylabel2use, 'FontSize', fSet.Fontsize_text)
set(gca,'YDir','reverse')
x_pos = get_value_range(tmp_x, 0.05);
y_pos = get_value_range(tmp_y, 0.15);
text(x_pos, y_pos, {betaString, pstring_chi, text_legend{1}}, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3])
% datatype for dose-response curve
colors = get_colors('spikewidth');
ncol = length(label_drug);
[fH, fSet] = plotj_initFig('width', 20, 'height', 14, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap.*[1.2 .8];
[P, h_text] = deal( NaN(length(label_drug), 2) ); % drugtype, unit-type,
[P_roc, h_text_roc] = deal( NaN(length(label_drug), 2) ); % drugtype, unit-type,
for idrug = 1:length(label_drug)
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
h_ax_subplot = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', idrug, 'axlabelDisplacement', [0.01, 0.02]);
title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
pos = get(h_ax_subplot, 'Position');
h_ax_inset = axes('Position', [pos(1)+0.28 pos(2)+0.035 0.1 0.11]) ; % inset
set(fH, 'currentaxes', h_ax_subplot);
unit2plot = idx_unit & att_table.unit_class==label_unitclass(iunitc);
idx_cond = att_table.attention=='Attend RF';
idx = unit2plot & idx_cond;
tmp_data = zeros(length(find(idx))/2,1);
tmp_data(:,1) = att_table.roc(idx & att_table.drug_on=='Drug off');
tmp_data(:,2) = att_table.roc(idx & att_table.drug_on=='Drug on');
% [bf10,p,CI,stats] = bf.ttest(tmp_data(:,1),tmp_data(:,2))
plotj_scatter(tmp_data, ...
'markerStyle', markerstyle, 'MarkerSize', markersize, ...
'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...
'markerEdgeColor', colors(iunitc,:), 'markerEdgeAlpha', 0.5, ...
P(idrug,iunitc) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');
p_string = get_significance_strings(P(idrug,iunitc), 'rounded', 0);
delta_data = tmp_data(:,2)-tmp_data(:,1);
d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');
fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
std(delta_data)/sqrt(size(tmp_data,1)), ...
x_pos = get_value_range(tmp_x, 0.05);
y_pos = get_value_range(tmp_y, 0.95 - 0.07 * (iunitc-1));
h_text(idrug,iunitc) = text(x_pos, y_pos, ...
sprintf('%s (n=%d)', p_string, size(tmp_data,1)), ...
'Color', colors(iunitc,:), 'FontSize', fSet.Fontsize_ax);
tmp_data(idx,:) = 1-tmp_data(idx,:);
% [bf10,p,CI,stats] = bf.ttest(tmp_data(:,1),tmp_data(:,2))
P_roc(idrug,iunitc) = compare_means(tmp_data(:,1),tmp_data(:,2), 1, 'rank');
bar_data = diff(tmp_data,1,2);
bar_y_err = std(bar_data)/sqrt(length(bar_data));
set(fH, 'currentaxes', h_ax_inset);
h_scatter = scatter(iunitc + (rand(length(bar_data),1)-0.5)*0.1, bar_data, ...
'.', 'MarkerFaceColor', colors(iunitc,:),'MarkerEdgeColor', colors(iunitc,:));
h_scatter.MarkerFaceAlpha = 0.5;
h_scatter.MarkerEdgeAlpha = 0.5;
plot(iunitc + [-0.25 0.25], [bar_y bar_y], 'Color', colors(iunitc,:), 'linew', fSet.LineWidth+1)
plot(iunitc + [0 0], bar_y + [-bar_y_err bar_y_err], 'Color', colors(iunitc,:), 'linew', fSet.LineWidth)
p_string = get_significance_strings(P_roc(idrug,iunitc), 'rounded', 0);
delta_data = tmp_data(:,2)-tmp_data(:,1);
d = computeCohen_d(tmp_data(:,2), tmp_data(:,1), 'paired');
fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
std(delta_data)/sqrt(size(tmp_data,1)), ...
% tmp_y = get(gca,'ylim');
%y_pos = get_value_range(tmp_y, 0.95);
h_text_roc(idrug,iunitc) = text(iunitc+0.1, y_pos, ...
'Color', colors(iunitc,:), ...
'FontSize', fSet.Fontsize_ax);
set(h_ax_inset, 'XTick', [1 2])
% ylabel(h_ax_inset, [plotj_symbol('Delta') ' | AUROC-0.5 |'],'FontSize', fSet.Fontsize_text)
ylabel(h_ax_inset, [plotj_symbol('Delta') ' AUROC-c'],'FontSize', fSet.Fontsize_text)
xlabel(h_ax_subplot, 'Attention AUROC no drug','FontSize', fSet.Fontsize_text)
ylabel(h_ax_subplot, 'Attention AUROC drug','FontSize', fSet.Fontsize_text)
end
AUROC: Dopamine, delta-AUROC 0.00 +- 0.02, p = 0.952, Cohens's d = 0.03
AUROC-CORRECT: Dopamine, delta-AUROC -0.01 +- 0.02, p = 0.542, Cohens's d = -0.16
AUROC: Dopamine, delta-AUROC -0.03 +- 0.01, p = 0.009, Cohens's d = -0.69
AUROC-CORRECT: Dopamine, delta-AUROC -0.02 +- 0.01, p = 0.080, Cohens's d = -0.38
AUROC: SCH23390, delta-AUROC -0.03 +- 0.04, p = 0.750, Cohens's d = -0.54
AUROC-CORRECT: SCH23390, delta-AUROC -0.03 +- 0.04, p = 0.750, Cohens's d = -0.54
AUROC: SCH23390, delta-AUROC -0.02 +- 0.04, p = 0.700, Cohens's d = -0.15
AUROC-CORRECT: SCH23390, delta-AUROC -0.04 +- 0.03, p = 0.278, Cohens's d = -0.33
[p_fdr, p_masked] = FDR(P, 0.05);
[p_fdr, p_masked_corr] = FDR(P_roc, 0.05);
plotj_text_emphasise(h_text, p_masked, 'italic');
plotj_text_emphasise(h_text, p_masked, 'bold');
% drugMI/attAUROC against ejection current
v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);
ylabel2use = 'Drug modulation index';
ylim2use = {[-0.5 0.5],[-0.4 0.2]};
idx_cond = att_table.attention=='Attend RF';
tmp_data = zeros(length(find(idx_cond))/2,1);
tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');
tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');
data2plot = diff(tmp_data,1,2);
ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];
ylim2use = {[-0.25 0.3],[-0.25 0.3]};
v_att_table = unstack(att_table, 'gain_log', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
% difference in attentional modulation of gain variabiliy
data2plot1 = (v_att_table.("Attend RF, Drug on") - v_att_table.("Attend away, Drug on")) ./ (v_att_table.("Attend RF, Drug on") + v_att_table.("Attend away, Drug on"));
data2plot2 = (v_att_table.("Attend RF, Drug off") - v_att_table.("Attend away, Drug off")) ./ (v_att_table.("Attend RF, Drug off") + v_att_table.("Attend away, Drug off"));
data2plot = data2plot1 - data2plot2;
idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});
ax(idrug) = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', iplot, 'axlabelDisplacement', [0.01, 0.02]);
title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
unit2plot = idx_unit & unitlist.unit_class==label_unitclass(iunitc);
xlabel_event = 'Ejection Current';
if length(unique(unitlist.unit_class(unit2plot)))==1
tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];
hscat = plotj_scatter(tmp_data, ...
'markerStyle', markerstyle, 'markerSize', markersize, ...
'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...
text_legend{1} = [label_unitclass{iunitc} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==iunitc))) ')'];
tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];
hscat = plotj_scatter(tmp_data, ...
'dataIndex', grp2idx(unitlist.unit_class(unit2plot)), ...
'markerStyle', markerstyle, 'MarkerSize', markersize, ...
'markerFaceColor', (colors(1:2,:)), 'markerFaceAlpha', [0.5 0.5], ...
'markerEdgeColor', (colors(1:2,:)), 'markerEdgeAlpha', [0.5 0.5], ...
text_legend{1} = [label_unitclass{1} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(1)))) ')'];
text_legend{2} = [label_unitclass{2} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(2)))) ')'];
data_table = table(unitlist.EjectCurrent_centered(unit2plot), data2plot(unit2plot), 'VariableNames', {'x', 'y'});
[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);
predict_data.x = predict_data.x + mean(unique(unitlist.EjectCurrent));
filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criterium));
filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criterium));
writetable(data_table, filename)
% fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)
if any(strcmp('fit_mean', predict_data.Properties.VariableNames))
h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);
cond1 = exist([filename(1:end-4) '_R.csv'], 'file');
cond2 = strcmpi(model,'quadratic');
utable = readtable([filename(1:end-4) '_R.csv']);
pstring_chi = [pstring_chi ' (U^{+})'];
pstring_chi = [pstring_chi ' (U^{-})'];
betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];
% model_color = colors(2,:);
betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];
xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)
ylabel(ylabel2use, 'FontSize', fSet.Fontsize_text)
% plot text indicating number of units
set(gca, 'Units', 'normalized');
x_pos = get_value_range(tmp_x, 0.6);
y_pos(1) = get_value_range(tmp_y, 0.8);
y_pos(2) = get_value_range(tmp_y, 0.9);
text(x_pos, y_pos(1), text_legend{1}, 'FontSize', fSet.Fontsize_text, 'Color', colors(1,:))
text(x_pos, y_pos(2), text_legend{2}, 'FontSize', fSet.Fontsize_text, 'Color', colors(2,:))
set(gca,'YDir','reverse')
x_pos = get_value_range(tmp_x, 0.05);
y_pos = get_value_range(tmp_y, 0.1);
text(x_pos, y_pos, {betaString, pstring_chi }, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3])
colors = get_colors('spikewidth');
selectivity_criteria = {'none', 'att', 'dru', 'att&dru'};
label_criteria = {'No subselection', 'Attention-selective', 'Drug-selective', {'Attention & Drug', 'selective'}};
ncol = length(selectivity_criteria);
nrow = length(label_drug);
[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap.*[1 .75];
idx_subplot = [1:ncol ; (ncol+1):ncol*2];
% drugMI/attAUROC against ejection current
v_drug_table = unstack(drug_table, 'MI', {'attention'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
data2plot = mean([v_drug_table.("Attend RF") v_drug_table.("Attend away")],2);
ylabel2use = 'Drug modulation index';
% data2plot = squeeze(diff(rate_ROC.roc,1,2));
% data2plot = squeeze(mean(data2plot(:,idx_attention_roc==1),2)); % average over relevant auroc conditions
idx_cond = att_table.attention=='Attend RF';
tmp_data = zeros(length(find(idx_cond))/2,1);
tmp_data(:,1) = att_table.roc(idx_cond & att_table.drug_on=='Drug off');
tmp_data(:,2) = att_table.roc(idx_cond & att_table.drug_on=='Drug on');
data2plot = diff(tmp_data,1,2);
ylabel2use = [plotj_symbol('Delta') ' attention AUROC'];
[P,h_text] = deal(zeros(length(selectivity_criteria), length(label_drug)));
for icrit = 1:length(selectivity_criteria)
selectivity_criterium = selectivity_criteria{icrit};
for idrug = 1:length(label_drug)
idx_unit = get_unit_selectivity(unitlist, selectivity_criterium, {'drug',label_drug(idrug)});
ax(idrug) = subtightplot(nrow, ncol, idx_subplot(iplot), fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', icrit);
% title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
h_title = title(label_criteria{icrit}, 'FontSize', fSet.Fontsize_title);
set(h_title, 'Units', 'normalized');
h_title.Position = h_title.Position+[0 0.2 0];
unit2plot = idx_unit & unitlist.unit_class==label_unitclass{iunitc};
xlabel_event = 'Ejection Current';
if length(unique(unitlist.unit_class(unit2plot)))==1
tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];
hscat = plotj_scatter(tmp_data, ...
'markerStyle', {'.'}, 'markerSize', markersize, ...
'markerFaceColor', colors(iunitc,:), 'markerFaceAlpha', 0.5, ...
text_legend{1} = [label_unitclass{iunitc} ' (n=' num2str(length(find(unit_class(idx_unit)==iunitc))) ')'];
tmp_data = [unitlist.EjectCurrent(unit2plot) data2plot(unit2plot)];
hscat = plotj_scatter(tmp_data, ...
'dataIndex', grp2idx(unitlist.unit_class(unit2plot)), ...
'markerStyle', {'o','o'}, 'MarkerSize', markersize, ...
'markerFaceColor', (colors(1:2,:)), 'markerFaceAlpha', [0.5 0.5], ...
'markerEdgeColor', (colors(1:2,:)), 'markerEdgeAlpha', [0.5 0.5], ...
text_legend{1} = [label_unitclass{1} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(1)))) ')'];
text_legend{2} = [label_unitclass{2} ' (n=' num2str(length(find(unitlist.unit_class(idx_unit)==label_unitclass(2)))) ')'];
data_table = table(unitlist.EjectCurrent(unit2plot), data2plot(unit2plot), 'VariableNames', {'x', 'y'});
[stats2report, predict_data, stats, model] = fitlme_singleVar_sequential(data_table, 0);
P(icrit,idrug) = stats2report.pValue;
filename = fullfile(path_population, sprintf('doseResponse_drugMI_%s_%s.csv', label_drug{idrug}, selectivity_criteria{icrit}));
filename = fullfile(path_population, sprintf('doseResponse_attAUROC_%s_%s.csv', label_drug{idrug}, selectivity_criteria{icrit}));
writetable(data_table, filename);
fprintf('Chi(%d) %1.2f, p = %1.3f\n', stats2report.deltaDF, stats2report.LRStat, stats2report.pValue)
if any(strcmp('fit_mean', predict_data.Properties.VariableNames))
h = plot(predict_data.x, predict_data.fit_mean, 'linewidth', fSet.LineWidth, 'color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_low, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
h = plot(predict_data.x, predict_data.fit_CI_high, '--', 'linewidth', fSet.LineWidth_in, 'Color', [0.3 0.3 0.3]);
[pstring_chi,starstring] = get_significance_strings(stats2report.pValue, 'rounded', 0);
cond1 = exist([filename(1:end-4) '_R.csv'], 'file');
cond2 = strcmpi(model,'quadratic');
utable = readtable([filename(1:end-4) '_R.csv']);
pstring_chi = [pstring_chi ' (U^{+})'];
pstring_chi = [pstring_chi ' (U^{-})'];
betaString = ['\beta' '_{1} = ' num2str(stats.Estimate, '%1.2e')];
% model_color = colors(2,:);
betaString = ['\beta' '_{2} = ' num2str(stats.Estimate, '%1.2e')];
if idrug==length(label_drug)
xlabel(xlabel_event, 'FontSize', fSet.Fontsize_text)
ylabel({label_drug{idrug},ylabel2use}, 'FontSize', fSet.Fontsize_text)
% plot text indicating number of units
set(gca, 'Units', 'normalized');
x_pos = get_value_range(tmp_x, 0.6);
y_pos(1) = get_value_range(tmp_y, 0.8);
y_pos(2) = get_value_range(tmp_y, 0.9);
text(x_pos, y_pos(1), text_legend{1}, 'FontSize', fSet.Fontsize_text, 'Color', colors(1,:))
text(x_pos, y_pos(2), text_legend{2}, 'FontSize', fSet.Fontsize_text, 'Color', colors(2,:))
set(gca,'YDir','reverse')
x_pos = get_value_range(tmp_x, 0.65);
y_pos = get_value_range(tmp_y, 0.05);
x_pos = get_value_range(tmp_x, 0.05);
y_pos = get_value_range(tmp_y, 0);
h_text(icrit,idrug) = text(x_pos, y_pos, {betaString, pstring_chi}, 'fontsize', fSet.Fontsize_text, 'color', [0.3 0.3 0.3]);
end
Chi(1) 0.42, p = 0.515
Chi(1) 0.82, p = 0.365
Chi(1) 4.71, p = 0.030
Chi(1) 0.74, p = 0.390
Chi(1) 1.03, p = 0.311
Chi(1) 0.80, p = 0.370
Chi(1) 0.69, p = 0.406
Chi(1) 0.80, p = 0.370
% multiple comparison correction, individual models
[p_fdr, p_masked] = FDR(P, 0.05);
plotj_text_emphasise(h_text, p_masked, 'italic', 2);
plotj_text_emphasise(h_text, p_masked, 'bold', 2);
[colors, labels] = get_colors('att_drug');
[colors_violin, labels_violin, labels_violin_short] = get_colors('att_drug');
colors_violin = reshape(colors_violin, [4,3]);
labels_violin = reshape(labels_violin, [4,1]);
ncol = length(label_drug);
[fH, fSet] = plotj_initFig('width',20,'height',15,'Journal',plot_conventions);
for idrug = 1:length(label_drug)
for iunitc = 1:length(unique(unitlist.unit_class))
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);
h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', idrug);
h_title = title([label_unitclass{iunitc} '-spiking'], 'FontSize', fSet.Fontsize_title);
gv = unstack(att_table(idx_unit,:), 'gv', {'attention_drug_on'}, 'GroupingVariables', {'unit','unit_class'}, 'VariableNamingRule', 'preserve');
variable_names = gv.Properties.VariableNames;
for iatt_drug = 1:length(labels_violin)
mus = [gv.(labels_violin{iatt_drug}).mus];
s2s = [gv.(labels_violin{iatt_drug}).s2s];
rhat = nanmean([gv.(labels_violin{iatt_drug}).rhat]);
gv_pop = GainVariability;
m2v = gv_pop.pred_variance(xval_gain);
h(iatt) = scatter(mus,s2s,fSet.MarkerSize*3,colors_violin(iatt_drug,:),'filled','MarkerFaceAlpha',0.5);
plot(m2v.mean,m2v.variance,'Color',colors_violin(iatt_drug,:),'linew',fSet.LineWidth);
xlim([0.8 max(xval_gain)])
ylim([0.8 max(xval_gain)])
plot([1e-10 max(s2s)*1.1],[1e-10 max(s2s)*1.1],'k','linew',1)
set(gca,'xscale','log','yscale','log')
% legend(h(:),labels(:),'Location','SouthEast','Fontsize',fSet.Fontsize_text)
if idrug==length(label_drug)
xlabel('Mean (spikes)','FontSize',fSet.Fontsize_text)
ylabel({[label_drug{idrug} label_drug_ext{idrug}],'Variance (spikes^{2})'},'FontSize',fSet.Fontsize_text)
pos = get(h_ax, 'Position');
h_ax_inset = axes('Position', [pos(1)+0.22 pos(2)+0.03 0.16 0.15]) ; % inset
idx_unit = get_unit_selectivity(att_table, selectivity_criterium, {'drug',label_drug(idrug)});
idx_unit = idx_unit & att_table.unit_class==label_unitclass(iunitc);
num_unit_plot = length(unique(att_table.unit(idx_unit)));
fprintf('\t Found %d units \n', num_unit_plot)
v_att_table = att_table(idx_unit, ismember(att_table.Properties.VariableNames, {'gain_log','drug_on','attention','unit','attention_drug_on'}));
v = violinplot(v_att_table.gain_log, v_att_table.attention_drug_on, 'ShowMean', true);
XTickLabel = get(gca,'XTickLabel');
v(iv).ScatterPlot.SizeData = 10;
idx = strcmpi(XTickLabel(iv), labels_violin);
v(iv).ViolinColor = colors_violin(idx,:);
h_ax_inset.XAxis.Visible = 'off';
% plot lines for each unit
[data_x, data_y] = deal(NaN(num_unit_plot, length(length(v))));
if ~any(isnan(v_att_table.gain_log))
% works if data is complete
data_x(:,iv) = v(iv).ScatterPlot.XData;
data_y(:,iv) = v(iv).ScatterPlot.YData;
% catch if nan values are present
idx = strcmpi(XTickLabel(iv), labels_violin);
idx_trial = v_att_table.attention_drug_on==labels_violin{idx};
data_y(:,iv) = v_att_table.gain_log(idx_trial);
plot(data_x', data_y', 'Color', [0.5 0.5 0.5 0.3]);
lme_population = fitlme(v_att_table,['gain_log ~ 1 + attention*drug_on + (1|unit)'],'DummyVarCoding','effects'); % fit interaction model with effect coding
disp(lme_population.Coefficients)
coefficients = lme_population.Coefficients;
betas = coefficients.Estimate(2:end);
betas_se = coefficients.SE(2:end);
p_val = coefficients.pValue(2:end);
% p_string = get_significance_strings(p_val, 'rounded', 0);
p_string = get_significance_strings(p_val, 'rounded', 0, 'factorstring', {'drug','att',sprintf('att%sdrug', '\times')});
p_string{end+1} = ['n = ' num2str(num_unit_plot)];
set(gca, 'ylim', YLIM .* [1 1.25])
x_pos = get_value_range(tmp_x, 0.95);
y_pos = get_value_range(tmp_y, 1.3);
h_text = text(x_pos, y_pos, p_string,'HorizontalAlignment','right');
ylabel('log(\sigma^{2}_G )', 'FontSize', fSet.Fontsize_text)
% [p_fdr, p_masked] = FDR(p_val, 0.05);
p_masked = [p_masked ; false];
plotj_text_emphasise(h_text, 1, 'italic', p_masked);
plotj_text_emphasise(h_text, 1, 'bold', p_masked);
end
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.6769 0.24109 -6.9553 47 9.5535e-09 -2.1619 -1.1919
{'drug_on_Drug on' } 0.23563 0.074521 3.1619 47 0.002744 0.085713 0.38555
{'attention_Attend RF' } -0.16931 0.074521 -2.272 47 0.027706 -0.31923 -0.019394
{'drug_on_Drug on:attention_Attend RF'} 0.097133 0.075694 1.2832 47 0.20571 -0.055145 0.24941
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.9154 0.1718 -11.149 63 1.5163e-16 -2.2587 -1.572
{'drug_on_Drug on' } 0.16573 0.045026 3.6809 63 0.00048392 0.075757 0.25571
{'attention_Attend RF' } -0.030159 0.045026 -0.66981 63 0.50543 -0.12014 0.059818
{'drug_on_Drug on:attention_Attend RF'} -0.031138 0.045026 -0.69156 63 0.49176 -0.12111 0.058839
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.1245 0.31895 -3.5257 8 0.0077812 -1.86 -0.38901
{'drug_on_Drug on' } 0.0032421 0.029576 0.10962 8 0.91541 -0.064961 0.071446
{'attention_Attend RF' } -0.1411 0.029576 -4.7706 8 0.0014075 -0.2093 -0.072895
{'drug_on_Drug on:attention_Attend RF'} 0.06162 0.029576 2.0834 8 0.070737 -0.0065831 0.12982
Fixed effect coefficients: DFMethod = 'Residual', Alpha = 0.05
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -1.6587 0.34668 -4.7845 37 2.7399e-05 -2.3611 -0.95625
{'drug_on_Drug on' } 0.31791 0.088241 3.6027 37 0.00092093 0.13911 0.4967
{'attention_Attend RF' } -0.1319 0.088241 -1.4948 37 0.14346 -0.31069 0.046895
{'drug_on_Drug on:attention_Attend RF'} 0.17933 0.087884 2.0405 37 0.048482 0.0012596 0.3574
mfactor = 1000; %multiplication_factor
colors = get_colors('att_drug');
nrow = length(label_drug);
[fH, fSet] = plotj_initFig('width', 10, 'height', 12, 'Journal',plot_conventions);
for idrug = 1:length(label_drug)
idx_rec = strcmpi(recordinglist.Drug, label_drug(idrug));
subtightplot(nrow, ncol, idrug, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
plotj_initAx(fSet, 'axlabel', idrug + 1, 'axlabelDisplacement', [0.07, 0.01]);
title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
[RT_table,RT_mat,colors_mat] = deal([]);
y = squeeze(mean(RT.RT(idx_rec,idrug_in,abs(idx_attention_cond-2)==iatt),3)) * mfactor;
text_legend{iattDrug} = sprintf('%s, %s', label_attention{iatt}, label_drug_onoff{idrug_in});
n = length(find(idx_rec));
drug = repmat(idrug_in,n,1);
RT_table = [RT_table ; table(rec,y,cond,drug,'VariableNames',{'recording','RT','attention','drug'})];
colors_mat = [colors_mat ; squeeze(colors(iatt,idrug_in,:))'];
RT_table.attention = categorical(RT_table.attention);
RT_table.drug = categorical(RT_table.drug);
v = violinplot(RT_mat, text_legend, 'ShowMean', true);
v(iv).ViolinColor = colors_mat(iv,:);
lme = fitlme(RT_table,['RT ~ 1 + (1|recording)'], 'DummyVarCoding', 'effects');% fit constant
lme2 = fitlme(RT_table,['RT ~ attention + (1|recording)'], 'DummyVarCoding', 'effects'); % fit linear model with attention condition
lme3 = fitlme(RT_table,['RT ~ attention + drug + (1|recording)'], 'DummyVarCoding', 'effects'); % fit linear model with states
lme4 = fitlme(RT_table,['RT ~ attention + drug + attention:drug + (1|recording)'], 'DummyVarCoding', 'effects'); % fit interaction
[BETA,BETANAMES,STATS] = fixedEffects(lme4);
comp1 = compare(lme, lme2, 'CheckNesting',true); % compare linear to constant
comp2 = compare(lme2, lme3, 'CheckNesting',true); %
comp3 = compare(lme3, lme4, 'CheckNesting',true); %
[p_fdr, p_masked] = FDR( [comp1.pValue(2) comp2.pValue(2) comp3.pValue(2)], 0.05);
pstring = get_significance_strings([STATS.pValue(2:end)], 'rounded', 0, 'factorstring', {'attention', 'drug', 'interaction'});
pstring{4} = ['n = ' num2str(length(find(idx_rec)))];
% plot text indicating number of units
set(gca, 'Units', 'normalized');
x_pos = get_value_range(tmp_x, 0.1);
y_pos = get_value_range(tmp_y, 0.85);
text(x_pos, y_pos, pstring, 'FontSize', fSet.Fontsize_text)
ylabel('Reaction time (ms)', 'FontSize', fSet.Fontsize_text)
if idrug == length(label_drug)
set(gca, 'xticklabelRotation', 25)
set(gca, 'xticklabel', '')
events = {'stim', 'cue', 'dim'};
xlabel_event = {'Time from stimulus onset (ms)', 'Time from cue onset (ms)', 'Time from first-dimming (ms)'};
xlim2use = {[-200 1000], [-200 1000], [-1000 200]};
ncol = length(events); % events
nrow = 2;%length(label_drug)-1;
[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap ./ [1 2];
for ievent = 1:length(events)
[pupil_timeseries, reclist] = get_population_data(recordinglist, 'pupil_timeseries', path_data, events{ievent});
if nrow==length(label_drug)
idx_rec = strcmpi(reclist.Drug, label_drug(idrug));
idx_rec = true(height(reclist),1);
nrecs = length(find(idx_rec));
pupil2plot = squeeze(nanmean(pupil_timeseries.pupil(idx_rec,:,:,:),3));
time = pupil_timeseries.time*1000;
timecenters = time(1):timestep:time(end);
p_drug = NaN(1,length(timecenters));
h_ax = subtightplot(nrow, ncol, iplot, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin);
plotj_initAx(fSet, 'axlabel', idrug, 'axlabelDisplacement', [0.07, 0.01]);
axh.YAxis.Visible = 'off'; % remove y-axis
y = squeeze(nanmean(pupil2plot(:,idrug_in,:),1));
y_err = squeeze(nanstd(pupil2plot(:,idrug_in,:),[],1))/sqrt(size(pupil2plot,1));
h(ilegend) = boundedline(...
'alpha', 'cmap', fSet.colors(idrug_in,:));
set(h(ilegend), 'LineWidth', fSet.LineWidth)
text_legend{ilegend} = sprintf('Drug %s', label_drug_onoff{idrug_in});
for itc = 1:length(timecenters)
timewin = timecenters(itc) + [-timewinsize/2 timewinsize/2];
idx_time = dsearchn(time(:), timewin(:));
tmp_data = squeeze(nanmean(pupil2plot(:,:,idx_time(1):idx_time(2)),3));
p_drug(itc) = compare_means(tmp_data(:,1), tmp_data(:,2), 1, 'rank');
[p_fdr, p_masked] = FDR( p_drug, 0.05);
plotj_brokenVector(timecenters, ylim2use(1)*0.95, p_masked, 'linewidth', 5, 'color', [0.5 0.5 0.5]);
hleg = legend(h, text_legend, 'autoupdate', 'off', 'FontSize', fSet.Fontsize_text, 'Box', 'Off');
plot([0 0], ylim, 'k','linew',1)
ylabel('Pupil diameter (z-score)', 'FontSize', fSet.Fontsize_text)
xlabel(xlabel_event{ievent}, 'FontSize', fSet.Fontsize_text)
title([label_drug{idrug} label_drug_ext{idrug}], 'FontSize', fSet.Fontsize_title)
text(h_ax, 500, -0.5, sprintf('n=%d',nrecs), 'FontSize', fSet.Fontsize_text)
end
Warning: Ignoring extra legend entries.
colors = get_colors('drug');
events = {'baseline', 'stim', 'cue', 'dim'};
event_title = {'Baseline','Stimulus evoked','Cue evoked','First-dimming'};
ncol = length(events); % events
[fH, fSet] = plotj_initFig('width', 20, 'height', 10, 'Journal',plot_conventions);
fSet.subplotGap = fSet.subplotGap/1.5;
[P, h_text] = deal(NaN(length(events), length(label_drug)));
for ievent = 1:length(events)
[pupil, reclist] = get_population_data(recordinglist, 'pupil_windows', path_data, events{ievent});
subtightplot(nrow, ncol, ievent, fSet.subplotGap, fSet.subplotMargin, fSet.subplotMargin)
plotj_initAx(fSet, 'axlabel', ievent+2, 'axlabelDisplacement', [0.05, 0.13]);
title(event_title{ievent}, 'FontSize', fSet.Fontsize_title)
text_p = cell(length(label_drug),1);
for idrug = 1:length(label_drug)
idx_rec = strcmpi(reclist.Drug, label_drug(idrug));
pupil2plot = squeeze(nanmean(pupil.pupil(idx_rec,:,:),3));
P(ievent,idrug) = compare_means(pupil2plot(:,1), pupil2plot(:,2), 1, 'rank');
text_p{idrug} = get_significance_strings(P(ievent,idrug), 'rounded', 0);
text_legend{idrug} = sprintf('%s (n=%d)', label_drug{idrug}, length(find(idx_rec)));
delta_data = pupil2plot(:,2)-pupil2plot(:,1);
d = computeCohen_d(pupil2plot(:,2), pupil2plot(:,1), 'paired');
fprintf('%s: %s, delta-%s %1.2f +- %1.2f, %s, Cohens''s d = %1.2f\n', ...
std(delta_data)/sqrt(length(find(idx_rec))), ...
pupil2plot = squeeze(nanmean(pupil.pupil,3));
h = plotj_scatter(pupil2plot, ...
'dataIndex', idx_rec+1, ...
'Markersize', markersize, 'markerStyle', {'o','s'}, ...
'MarkerEdgeColor', colors, 'MarkerFaceColor', colors, ...
'MarkerEdgeAlpha', [0.5 0.5] ,'MarkerFaceAlpha', [0.5 0.5]);
% set(gca,'ydir','reverse')
x_pos = get_value_range(tmp_x, 0.6);
for idrug = 1:length(label_drug)
y_pos = get_value_range(tmp_y, 0.2-0.1*(idrug-1));
h_text(ievent,idrug) = text(x_pos, y_pos, text_p{idrug}, 'FontSize', fSet.Fontsize_text, 'color', colors(idrug,:));
ylabel({'Drug (z-score)'}, 'FontSize', fSet.Fontsize_text)
xlabel({'No drug (z-score)'}, 'FontSize', fSet.Fontsize_text)
hleg = legend(h, text_legend, 'FontSize', fSet.Fontsize_text, 'Box', 'Off', 'Location', 'NorthEast');
hleg.Position = hleg.Position + [0.05 0.3 0 0];
end
PUPIL: Dopamine, delta-pupil 0.02 +- 0.04, p = 0.821, Cohens's d = 0.13
PUPIL: SCH23390, delta-pupil 0.01 +- 0.05, p = 0.890, Cohens's d = 0.04
PUPIL: Dopamine, delta-pupil 0.10 +- 0.02, p < 0.001, Cohens's d = 1.09
PUPIL: SCH23390, delta-pupil 0.10 +- 0.03, p = 0.004, Cohens's d = 0.79
PUPIL: Dopamine, delta-pupil 0.08 +- 0.03, p = 0.051, Cohens's d = 0.66
PUPIL: SCH23390, delta-pupil 0.07 +- 0.03, p = 0.031, Cohens's d = 0.62
PUPIL: Dopamine, delta-pupil 0.05 +- 0.03, p = 0.193, Cohens's d = 0.37
PUPIL: SCH23390, delta-pupil 0.06 +- 0.03, p = 0.023, Cohens's d = 0.54
Warning: Ignoring extra legend entries.
[p_fdr, p_masked] = FDR(P, 0.05);
plotj_text_emphasise(h_text, p_masked, 'italic');
plotj_text_emphasise(h_text, p_masked, 'bold');